System and method of modeling erythropoiesis including iron homeostasis

ABSTRACT

A method of adjusting a patient&#39;s undesired hematocrit and/or hemoglobin concentration to a value within a desired range at a predetermined time with an erythropoiesis stimulating agent (ESA) regimen includes obtaining patient parameters required for input into a model for predicting the patient&#39;s hematocrit and/or hemoglobin concentration at a predetermined time with a selected ESA administration regimen, the model including iron homeostasis, and employing the patient parameters and an initially selected ESA administration regimen in the model to predict the patient&#39;s hematocrit and/or hemoglobin concentration at the predetermined time with the initially selected ESA administration regimen. The method then includes administering ESA to the patient with an ESA administration regimen predicted to adjust the patient&#39;s hematocrit and/or hemoglobin concentration to the desired range at the predetermined time. The method can be implemented in a hemodialysis system.

RELATED APPLICATION(S)

This application claims the benefit of U.S. Provisional Application No.61/722,503, filed on Nov. 5, 2012. The entire teachings of thisapplication are incorporated herein by reference.

BACKGROUND OF THE INVENTION

Red blood cells (erythrocytes) are essential for the transport of oxygenthrough the body. An understanding of the regulation of red blood cellproduction, called erythropoiesis, is important for the treatment ofpatients in a variety of clinical situations. Patients that arescheduled for elective surgery, such as hip or transplant surgery, canbe prescribed an erythropoiesis stimulating agent (ESA) to compensatefor the expected loss of blood, thus obviating the need for allogenicblood transfusions by raising the patient's hematocrit and/or hemoglobinconcentration to a desired range at the predetermined time, inexpectation of the surgery. Erythropoiesis stimulating agents (ESAs),including recombinant human erythropoietin, exert hematological effectsanalogous to the hormone erythropoietin (EPO), which is released intothe blood stream by the kidneys based on a negative feedback mechanismthat reacts to the partial pressure of oxygen in the blood. ESAtreatment regimens are also prescribed for patients who suffer frominsufficient erythropoiesis, such as cancer patients recovering from theeffects of chemotherapy, and chronic kidney disease patients whosekidneys can no longer produce sufficient amounts of EPO. The dose andfrequency of administration of an ESA treatment regimen are oftendetermined based on the prior experience of the physician and onestablished guidelines, because predictive models of erythropoiesisunder an ESA treatment regimen are not readily available.

Therefore, there is a need for a reliable predictive model oferythropoiesis under various ESA treatment regimens.

SUMMARY OF THE INVENTION

A model for erythropoiesis in humans is provided that is based onstructured population models for the different cell stages indevelopment, from stem cells in bone marrow to erythrocytes in the bloodstream.

In one embodiment, a method of adjusting a patient's hematocrit and/orhemoglobin concentration to a desired range at a predetermined time withan erythropoiesis stimulating agent (ESA) regimen includes obtainingpatient parameters required for input into a model for predicting thepatient's hematocrit and/or hemoglobin concentration at a predeterminedtime with a selected ESA administration regimen, the model includingiron homeostasis, and employing the patient parameters and an initiallyselected ESA administration regimen in the model to predict thepatient's hematocrit and/or hemoglobin concentration at thepredetermined time with the initially selected ESA administrationregimen. Optionally, if the patient's hematocrit and/or hemoglobinconcentration is not predicted by the model to be in the desired rangeat the predetermined time, the method includes employing the model withone or more different ESA administration regimens until the modelpredicts that the patient's hematocrit and/or hemoglobin concentrationwill be in the desired range at the predetermined time. The method thenincludes administering ESA to the patient with the ESA administrationregimen predicted to adjust the patient's hematocrit and/or hemoglobinconcentration to the desired range at the predetermined time. The ESAadministration regimen can include iron administration, and the modelincludes iron homeostasis with iron administration, and modeling theinfluence of iron homeostasis on erythropoiesis. The patient parameterscan include the starting hemoglobin concentration in the patient'sblood, the total blood volume of the patient, the lifespan of red bloodcells (RBCs) of the patient, the mean corpuscular volume of the RBCs,the rate of neocytolysis in the patient's blood, the distribution ofcell hemoglobin content in the patient's blood, the distribution of redblood cell size, iron storage level of the patient, C-reactive protein(CRP) levels, interleukin-6 (IL-6) levels, and the hepcidin levels ofthe patient. The predetermined time can be, for example, in a range ofbetween about 5 days and about 200 days into the ESA administrationregimen. In some embodiments, the patient undergoes a medical procedureprior, during, or after initiation of an ESA administration regimen,including medical procedures such as blood donation, surgery, anddialysis, or any combination thereof. For dialysis patients, the desiredhematocrit can be, for example, in the range of between about 28 percentand about 36 percent and the desired hemoglobin concentration can be,for example, in a range of between about 9.5 g/dL and about 12 g/dL.

In yet another embodiment, a computer system for adjusting a patient'shematocrit and/or hemoglobin concentration to a desired range at apredetermined time with an erythropoiesis stimulating agent (ESA)regimen includes a user input means for determining patient parametersfrom a user, a digital processor coupled to receive determined patientdata from the input means, wherein the digital processor executes amodeling system in working memory, and an output means coupled to thedigital processor, the output means provides to the user the patient'shematocrit and/or hemoglobin concentration under each ESA administrationregimen at the predetermined time. The modeling system employs thepatient parameters and an initially selected ESA administration regimenin the model to predict the patient's hematocrit and/or hemoglobinconcentration at the predetermined time with the initially selected ESAadministration regimen, the model including iron homeostasis.Optionally, if the patient's hematocrit and/or hemoglobin concentrationis not predicted by the model to be in the desired range at thepredetermined time, the model can be employed with one or more differentESA administration regimens until the model predicts that the patient'shematocrit and/or hemoglobin concentration will be in the desired rangeat the predetermined time.

In still another embodiment, a method of determining an ESAadministration regimen to bring a patient's hematocrit and/or hemoglobinconcentration into a desired range at a predetermined time with anerythropoiesis stimulating agent (ESA) regimen includes obtainingpatient parameters required for input into a model for predicting thepatient's hematocrit and/or hemoglobin concentration at a predeterminedtime with a selected ESA administration regimen, the model includingiron homeostasis, and employing patient parameters and an initiallyselected ESA administration regimen in the model to predict thepatient's hematocrit and/or hemoglobin concentration at thepredetermined time with the initially selected ESA administrationregimen. Optionally, if the patient's hematocrit and/or hemoglobinconcentration is not predicted by the model to be in the desired rangeat the predetermined time, the method includes employing the model withone or more different ESA administration regimens until the modelpredicts that the patient's hematocrit and/or hemoglobin concentrationwill be in the desired range at the predetermined time. The method thencan include administering ESA to the patient with the ESA administrationregimen predicted to adjust the patient's hematocrit and/or hemoglobinconcentration to the desired range at the predetermined time.

In another embodiment, a method of confirming that an ESA administrationregimen is maintaining a patient's hematocrit and/or hemoglobinconcentration within a desired range at a predetermined time with anerythropoiesis stimulating agent (ESA) regimen includes obtainingpatient parameters required for input into a model for predicting thepatient's hematocrit and/or hemoglobin concentration at a predeterminedtime with a selected ESA administration regimen, the model includingiron homeostasis, and employing patient parameters and an initiallyselected ESA administration regimen in the model to predict thepatient's hematocrit and/or hemoglobin concentration at thepredetermined time with the initially selected ESA administrationregimen. Optionally, if the patient's hematocrit and/or hemoglobinconcentration is not predicted by the model to be in the desired rangeat the predetermined time, the method includes employing the model withone or more different ESA administration regimens until the modelpredicts that the patient's hematocrit and/or hemoglobin concentrationwill be in the desired range at the predetermined time. The method thencan include administering ESA to the patient with the ESA administrationregimen predicted to adjust the patient's hematocrit and/or hemoglobinconcentration to the desired range at the predetermined time.

In yet another embodiment, a hemodialysis system includes a hemodialysismachine including a blood pump, a modular drug delivery devicecomprising at least one erythropoiesis stimulating agent (ESA) pump, atleast one iron replacement product pump, at least one ESA vial holder,and at least one iron replacement product vial holder, and a housing towhich the blood pump and the modular drug delivery device are secured.

The hemodialysis system further includes a blood line set comprisingmultiple blood lines and a drip chamber in fluid communication with themultiple blood lines, at least one of the blood lines being connected tothe drip chamber, the blood line set being connected to the blood pumpin a manner such that, when the blood line set is connected to a patientand the blood pump is operated, blood of the patient is passed throughthe blood line set. Additionally, the hemodialysis system includes adrug administration fluid line set comprising multiple drug lines, atleast two of the drug lines of the drug line set being connected to thedrip chamber of the blood line set, and the drug line set beingconnected to the at least one ESA pump and the at least one ironreplacement product pump in a manner such that, when the drug line setis fluidly connected to an ESA vial contained in the at least one ESAvial holder and to an iron replacement product vial contained in the atleast one iron replacement product vial holder, and the at least one ESApump and the at least one replacement product pump are operated, ESA andiron replacement product are delivered from the ESA vial and the ironreplacement vial, respectively, to the drip chamber of the blood lineset via the drug line set.

The hemodialysis system is controlled by a control unit configured toadjust the patient's undesired hematocrit and/or hemoglobinconcentration to a value within a desired range at a predetermined timeby i) employing a model for predicting the patient's hematocrit and/orhemoglobin concentration at a predetermined time, the model includingiron homeostasis, and by ii) operating the blood pump, the at least oneESA pump, and the at least one iron replacement product pump to deliverblood, ESA, and iron replacement product to the drip chamber via the atleast one of the blood lines connected to the drip chamber and the atleast two of the drug lines connected to the drip chamber, according toan ESA administration regimen predicted to adjust the patient'shematocrit and/or hemoglobin concentration to a value within the desiredrange at the predetermined time.

This invention has many advantages, including the achievement of an ESAregimen needed for a desired hematocrit and/or hemoglobin concentrationfor a patient, thereby, on the one hand, alleviating insufficienterythropoiesis, and, on the other hand, preventing excessively high ESAdose levels that raise the patient's blood pressure and increase thepatient's risk of stroke and cardiovascular disease.

BRIEF DESCRIPTION OF THE DRAWINGS

The foregoing will be apparent from the following more particulardescription of example embodiments of the invention, as illustrated inthe accompanying drawings in which like reference characters refer tothe same parts throughout the different views. The drawings are notnecessarily to scale, emphasis instead being placed upon illustratingembodiments of the present invention.

FIG. 1 is a schematic illustration of the influence of EPO and iron oncell stages during erythropoiesis in a human being.

FIG. 2 is an illustration of a hemodialysis machine that includes amodular drug delivery device. A drug administration fluid line set andmultiple drug vials are secured to the modular drug delivery device.

FIG. 3 is a perspective view of the modular drug delivery device of thehemodialysis machine shown in FIG. 2 and the drug administration fluidline set and drug vials that are secured to the modular drug deliverydevice.

FIG. 4 is a schematic illustration of the iron cycle in an adult. Thenumbers indicate the amount of iron (in mg) that is shifted per dayunder normal conditions. Approximately 1 mg of iron is excreted dailybecause of shedding of epithelial cells, via urine and sweat. Womenadditionally lose iron due to menses or pregnancy.

FIG. 5 is a schematic illustration of systemic iron regulation.

FIG. 6A is a schematic illustration of the organizational diagram of themodel. The feedback for the iron submodel is not illustrated here, forreasons of clarity.

FIG. 6B is a schematic illustration of the organizational diagram of themodel including the feedback for the iron submodel.

FIG. 7 is a schematic illustration of the fluxes through the differentboundaries. Parts of the boundaries which are non-zero are depicted ingray.

FIG. 8 is a schematic illustration of the organizational diagram of theiron submodel.

FIG. 9 is a list of definitions for model parameters for the equationsset forth in Section 4.

FIGS. 10-1-6 are graphs of model results of the total amount [mg] ofiron as a function of time [days] of FIG. 10-1: Plasma, FIG. 10-2:Erythrocytes, FIG. 10-3: Storage, FIG. 10-4: Precursor Cells, FIG. 10-5:Macrophages, FIG. 10-6: Total iron for the iron model for a 70 kg manstarted near the steady state.

FIGS. 11-1-6 are graphs of model results of the total amount [mg] ofiron as a function of time [days] of FIG. 11-1: Plasma, FIG. 11-2:Erythrocytes, FIG. 11-3: Storage, FIG. 11-4: Precursor Cells, FIG. 11-5:Macrophages, FIG. 11-6: Total Iron for the iron model for a 70 kg manstarted with increased amounts of iron incorporated in hemoglobin anddecreased storage.

FIG. 12 is a graph of transferrin receptors (#×10⁴) as a function ofcell age (days).

FIGS. 13 and 14A-14B are graphs of hemoglobin as a function of thecalendric age of precursor cell, i.e., days since cell became aprecursor cell. Depicted are a period of 6, 7 and 8 days.

FIG. 15 is a graph of hemoglobin (pg) as a function of days since cellbecame a precursor cell (days), for a period of 7 days. Shown are areference curve (solid line), that reflects the average hemoglobinaccumulation in cells and two curves (dashed lines) that depicttolerable variations in this average, thereby forming a corridor ofacceptable hemoglobin values around the average which acts as areference curve.

FIG. 16 is a graph of the rate of apoptosis as a function of hemoglobinconcentration (g/dl) for a certain cell age.

FIG. 17 illustrates the family of functions

(h,η). It is a graph of the rates of apoptosis as a function ofhemoglobin concentration (g/dl) and calendric age (days). Forexemplification, only calendric ages between 4-5 days are depicted. Thesolid line illustrates the hemoglobin reference curve, and thedash-dotted lines illustrate the corridor around this line (see FIG.15).

FIG. 18 is a graph describing the rates of apoptosis as a function ofhemoglobin concentration (g/dl) and calendric age illustrating theevaluation of the functions

(h,η) at the grid points η_(i) and the corresponding hemoglobin valuesh_(i). The resulting values (large dots) are interpolated using cubicsplines. For exemplification, only calendric ages between 4-5 days aredepicted.

FIG. 19 is a graph of the rate of apoptosis as a function of calendricage (days), illustrating the computation of the apoptosis function α.The values of the functions

at the discrete time points η_(i) and the corresponding hemoglobinvalues h_(i) (large dots) are interpolated using cubic splines (dashedline). The resulting curve represents the rate of apoptosis of thecells. The function is only plotted for a calendric age in a range ofbetween 4 days and 5 days, for clarity.

FIG. 20 is a schematic view of a computer network in which the presentinvention can be implemented.

FIG. 21 is a block diagram of a computer of the network of FIG. 20.

DETAILED DESCRIPTION OF THE INVENTION

Red blood cells (erythrocytes) are essential for the distribution ofoxygen through the body to organs and tissues. They take up oxygen inthe lungs and deliver it to tissues while squeezing through thecapillaries. To fulfill this task properly, they are highly specialized.For instance, being shaped like biconcave disks optimizes the oxygenexchange. Furthermore, the nuclei, organelles, and mitochondria havebeen expelled in earlier developmental stages, and therefore more spacefor hemoglobin, the molecule which oxygen binds to, is available.Erythrocytes are very deformable and can therefore pass capillaries halftheir diameter. During microcirculation, they have to withstand highshear stresses, rapid elongation, folding, and deformation. Over time,the cell membrane is damaged by these extraordinary stresses. Because ofthe lack of nuclei and organelles, red blood cells cannot divide orrepair their cell membranes. Senescent erythrocytes lose theirflexibility due to their fragmented membranes. These stiff cells coulddo harm to small capillaries or even clog them. To avoid this potentialharm, old erythrocytes are recognized by phagocytes and destroyed. Thisphagocytosis mainly takes place in the spleen and cells of thereticulo-endothelial system (RES). See Jandl, J. H., Blood: Textbook ofHematology, 2^(nd) Ed. Little, Brown and Company, 1996 (hereinafter“Jandl”).

To compensate for phagocytosis of senescent red blood cells, it isnecessary to build new erythrocytes continuously. The maturation ofundifferentiated stem cells to mature erythrocytes is callederythropoiesis and takes place in the bone marrow. Erythropoiesis notonly has to compensate for the continuous loss of old erythrocytes, butalso for the additional loss of cells due to random breakdown, as wellas due to internal and external bleeding. Furthermore, the number of redblood cells has to be adjusted to varying environmental conditions, asfor instance a transition from low to high altitudes or vice versa, byincreasing the rate of erythropoiesis, or, conversely, by neocytolysis,a process believed to be wherein macrophages start to phagocytose youngerythrocytes (neocytes).

During the process of erythropoiesis, the cell population undergoes aseries of proliferations and differentiations. Starting frommultipotential stem cells, erythroid cells mature to BFU-Es (earlieststage of erythroid committed cells), CFU-Es, different stages oferythroblasts, and finally reticulocytes. The reticulocytes are releasedfrom the bone marrow into blood and mature within 1-2 days toerythrocytes (see FIG. 1).

The primary control of erythropoiesis is governed by the hormoneerythropoietin (EPO). EPO is released into the blood stream by thekidneys based on a negative feedback mechanism that reacts to thepartial pressure of oxygen in blood. The concentration of EPO affectsthe number of circulating red blood cells by determining the number ofcells that mature into erythrocytes, either by recruitment or bypreventing apoptosis (programmed cell death), and by affecting thevelocity of maturing of progenitor and precursor cells. Thus,disturbances in oxygen delivery can be adjusted for by an adaptiveresetting of the rate of erythropoiesis. Additionally, as alreadymentioned above, there exists a physiological process which affects theselective degradation of young erythrocytes in situations of red cellexcess, called neocytolysis. Neocytolysis seems to be triggered by adrop in the EPO level. See Rice, L., and Alfrey, C. P., CellularPhysiology and Biochemistry, 2005, Vol. 15, pp. 245-250 (hereinafter“Rice 2005”); Rice, L., Alfrey, C. P., Driscoll, T., Whitley, C. E.,Hachey, D. L., and Suki, W., Amer. J. Kidney Diseases, 1999, Vol. 33,pp. 59-62 (hereinafter “Rice 1999”); and Rice, L., W. Ruiz, W.,Driscoll, T., Whitley, C. E., Tapia, R., Hachey, D. L., Conzales, G. F.,and Alfrey, C. P., Annals Internal Medicine. 2001, Vol. 134, pp. 652-656(hereinafter “Rice 2001”).

Another critical factor for effective erythropoiesis is the availabilityof iron which is indispensable for hemoglobin synthesis. If the body isnot able to provide sufficient iron for this process, then ineffectiveerythropoiesis will result. See Finch, S., Haskins, D., and Finch, C.A., The Journal of Clinical Investigation. 1950, Vol. 29, pp. 1078-1086;and Lichtman, M. A., E. Beutler, E., Kipps, T. J., Seligsohn, U.,Kaushansky, K., and Prchal, J. T. (editors), Williams Hematology, 7thedition, New York, McGraw-Hill, 2005 (hereinafter “WilliamsHematology”). In normal subjects, the total iron content of the bodystays within narrow limits (iron overload is toxic). Once an atom ofiron enters the body it is conserved with remarkable efficiency and canremain in the body for more than ten years. Iron is lost via loss ofcells (especially epithelial cells), bleeding and loss of very smallamounts via urine and sweat. The balance of iron content is achieved byabsorption and not by control of excretion. If the plasma concentrationof iron is too low, then the level of the hormone hepcidin is decreased.The consequence of a lower hepcidin level is that more iron is taken upvia the duodenum and more iron is released from macrophages and from thestores. See Crichton, R., Iron Metabolism. From Molecular Mechanisms toClinical Consequences, New York, J. Wiley, 2009 (hereinafter“Crichton”); Fleming, M. D., J. Amer. Soc. Hematology, 2008, pp. 151-158(hereinafter “Fleming”). Patients suffering from inflammation, such asdialysis patients, typically have higher hepcidin levels. Increasingiron availability in inflamed dialysis patients can be achieved by anincrease of parenteral iron by increasing dose, frequency, or both, andby reducing inflammation by diagnosis and treatment of sources ofinflammation, e.g., barrier breakdown (i.e., skin, periodontal disease,intestinal congestion), pulmonary or urinary tract infection, thrombosedfistulas or catheter, and by subsequent specific therapy, e.g.,antibiotics, catheter removal, aseptic techniques when manipulatingin-dwelling catheters, and surgical debridement of skin ulcers.

Since individual cells in the various cell populations which have to beconsidered have to be distinguished according to their age,age-structured population models are needed in order to describe thedevelopment of the cell populations. Besides these age-structuredpopulation models, the model of this invention includes a feedback loopincluding erythropoietin.

In one embodiment, a method of adjusting a patient's hematocrit and/orhemoglobin concentration to a desired range at a predetermined time withan erythropoiesis stimulating agent (ESA) regimen includes obtainingpatient parameters required for input into a model for predicting thepatient's hematocrit and/or hemoglobin concentration at a predeterminedtime with a selected ESA administration regimen, the model includingiron homeostasis, and employing the patient parameters and an initiallyselected ESA administration regimen in the model to predict thepatient's hematocrit and/or hemoglobin concentration at thepredetermined time with the initially selected ESA administrationregimen. Examples of ESAs are provided in Table 1 (adapted fromPhurrough S, Jacques L, Ciccanti M, Turner T, Koller E, Feinglass S:“Proposed Coverage Decision Memorandum for the Use of ErythropoiesisStimulating Agents in Cancer and Related Neoplastic Conditions”; Centersfor Medicare and Medicaid Services; Administrative File: CAG #000383N;May 14, 2007). See also Pfeffer, M. A., Burdmann, E. A., Chen, C-Y.,Cooper, M. E., de Zeeuw, D., Eckardt, K-U., Feyzi, J. M., Ivanovich, P.,Kewalramani, R., Levey, A. S., Lewis, E. F., McGill, J. B., McMurray, J.J. V., Parfrey, P., Parving, H-H., Remuzzi, G., Singh, A. K., Solomon,S. D., Toto, R., The New England Journal of Medicine, 2009, 361(21), pp.2019-2032; and Singh, A. K., Szczech, L., Tang, K. L., Barnhart, H.,Sapp, S., Wolfson, M., Reddan, D., The New England Journal of Medicine,2006, 355(20), pp. 2085-2098 (hereinafter “Singh et al.”). Note thatunlike other ESAs listed in Table 1, Peginesatide is not a biologicallyderived EPO; it is a synthetic peptide that stimulates EPO receptors.

TABLE 1 Erythropoiesis Stimulating Agents: EPO = erythropoietin.Compound Drug Names Manufacturer EPO α Epogen ® Amgen EPO α Procrit ®Amgen EPO α (w/o serum Eprex ® Epypo ® J&J subsidiary (Otho albumin)Epopen ® Epoxitin ® Biologics) Globuren ® EPO β (Neo)Recormon Roche EPOβ Erantin ® Boehringer Mannheim (Spain), Roche (Spain) EPO β Epoch ®Chugai EPO δ in human cell Dynepo Gene Activated Aventis Transkaryoticlines EPO Therapies EPO Ω Epomax ® Hemax ® Baxter Hemax ®-Eritron ®Modified EPO α Aranesp ® Amgen Darbepoietin Modified EPO α Nespo ® AmgenDarbepoietin Modified EPO β Mircera ® Roche Continuous EPO ReceptorActivator (Pegylation) Peginesatide Omontys ® Affymax

Optionally, if the patient's hematocrit and/or hemoglobin concentrationis not predicted by the model to be in the desired range at thepredetermined time with the initially selected ESA administrationregimen, the method includes employing the model with one or moredifferent ESA administration regimens until the model predicts that thepatient's hematocrit and/or hemoglobin concentration will be in thedesired range at the predetermined time. The method then includesadministering ESA to the patient with an ESA administration regimenpredicted to adjust the patient's hematocrit and/or hemoglobinconcentration to a value within the desired range at the predeterminedtime. The ESA administration regimen can include iron administration,and the model includes iron homeostasis with iron administration, andmodeling the influence of iron homeostasis on erythropoiesis. Thepatient parameters can include, for example, the starting hematocritand/or hemoglobin concentration in the patient's blood, the total bloodvolume of the patient, the lifespan of red blood cells (RBCs) of thepatient, the mean corpuscular volume of the RBCs, the rate ofneocytolysis in the patient's blood, the distribution of cell hemoglobincontent in the patient's blood, the distribution of red blood cell size,iron storage level of the patient, transferrin saturation (TSAT), serumtransferrin receptors, red cell zinc protoporphyrin level, hypochromicred blood cells, hypochromic reticulocytes, ferritin, endogenous EPOlevels, C-reactive protein (CRP) levels, interleukin-6 (IL-6) levels,and the hepcidin levels of the patient.

The starting hematocrit and/or hemoglobin concentration in the patient'sblood can be obtained from routine laboratory measurements known in theart. The total blood volume (BV) of the patient can be estimated asdescribed further below, or measured by use of radio-labeling red bloodcells with chromium-51 to estimate red blood cell volume (RCV) and usingthe formulaBV=RCV/(0.9*Hctv)where Hctv is the venous hematocrit, obtained from routine laboratorymeasurements known in the art. See Albert S. N., Blood volumemeasurement, In: Nuclear Medicine In Vitro. 2 ed. Philadelphia: JBLippincott Co., 1983; Bernard P. J., Nouv Rev Fr Hematol 1994, 36(2),pp. 155-157; and International Committee for Standardization inHaematology: Recommended methods for measurement of red-cell and plasmavolume, J Nucl Med 1980, 21(8), pp. 793-800. The lifespan of RBCs of thepatient can be estimated from endogenous alveolar carbon monoxideconcentrations. See Strocchi A., Schwartz S., Ellefson M., Engel R. R.,Medina A., Levitt M. D., J Lab Clin Med 1992, 120(3), pp. 392-399. Themean corpuscular volume can be obtained from routine laboratorymeasurements known in the art. The rate of neocytolysis in the patient'sblood can be estimated from correlations with reduced expression of CD44(homing-associated cell adhesion molecule) and CD71 (transferrinreceptor). See Chang C. C., Chen Y., Modi K., Awar O., Alfrey C., RiceL., J Investig Med 2009, 57(5), pp. 650-654.

The models of this invention can track the patient's predictedhematocrit and/or hemoglobin concentration over time, such as betweenabout 5 days and about 200 days of the ESA administration regimen. Thepredetermined time can be any future time after an ESA administrationregimen is selected and the predicted regimen is initiated. In someembodiments, the patient undergoes a medical procedure prior, during, orafter the ESA administration regimen, such as blood donation, surgery,and dialysis, or any combination thereof. For dialysis patients, thedesired hematocrit is typically in the range of between about 28 percentand about 36 percent and the desired hemoglobin concentration istypically in a range of between about 9.5 g/dL and about 12 g/dL. SeeDrüeke, T. B., Locatelli, F., Clyne, N., Eckardt, K-U., Macdougall, I.C., Tsakiris, D., Burger, H-U., Scherhad, A., The New England Journal ofMedicine, 2006, 355(20), pp. 2071-2084; Parfrey, P. S., 2006, AmericanJournal of Kidney Diseases, 47(1), pp. 171-173; Strippoli, G. F. M.,Craig, J. C., Manno, C., Schena, F. P., Journal of the American Societyof Nephrology, 2004, 15, pp. 3154-3165 (hereinafter “Strippoli et al.”);Volkova, N., Arab, L., Evidence-Based Systematic Literature Review ofHemoglobin/Hematocrit and All-Cause Mortality in Dialysis Patients,2006, 47(1), pp. 24-36. For elective orthopaedic surgery patients, thedesired hemoglobin concentration for males and females is typicallygreater than or equal to 13 g/dL, and 12 g/dL, respectively. SeeGoodnough L. T., Maniatis A., Earnshaw P., Benoni G., Beris P., BisbeE., Fergusson D. A., Gombotz H., Habler O., Monk T. G., Ozier Y,Slappendel R., and Szpalski M., British Journal of Anaesthesia 106 (1)pp. 13-22 (2011).

In another embodiment, a method of determining a patient's hematocritand/or hemoglobin concentration within a desired range at apredetermined time with an erythropoiesis stimulating agent (ESA)regimen includes obtaining patient parameters required for input into amodel for predicting the patient's hematocrit and/or hemoglobinconcentration at a predetermined time with a selected ESA administrationregimen, the model including iron homeostasis, and employing the patientparameters and an initially selected EPO administration regimen in themodel to predict the patient's hematocrit and/or hemoglobinconcentration at the predetermined time with the initially selected ESAadministration regimen. Optionally, if the patient's hematocrit and/orhemoglobin concentration is not predicted by the model to be in thedesired range at the predetermined time, or a different ESAadministration regimen is desired due to other considerations, themethod includes employing the model with one or more different ESAadministration regimens until the model predicts that the patient'shematocrit and/or hemoglobin concentration will be in the desired rangeat the predetermined time. The method then can include administering ESAto the patient with the ESA administration regimen predicted to adjustthe patient's hematocrit and/or hemoglobin concentration to a valuewithin the desired range at the predetermined time.

In yet another embodiment, a hemodialysis system includes a hemodialysismachine including a blood pump, a modular drug delivery devicecomprising at least one erythropoiesis stimulating agent (ESA) pump, atleast one iron replacement product pump, at least one ESA vial holder,and at least one iron replacement product vial holder, and a housing towhich the blood pump and the modular drug delivery device are secured.

Referring to U.S. Pat. No. 8,382,696 B2 issued to Beiriger et al. onFeb. 26, 2013, and as shown in FIG. 2, a hemodialysis system 100includes a hemodialysis machine 101 that has a drug delivery system 102.The drug delivery system 102 includes a modular drug delivery device 103and a disposable drug administration fluid line set 109 that isconnected to the drug delivery device 103. A drug delivery line 104 ofthe drug administration fluid line set 109 is fluidly connected to ablood circuit of the hemodialysis system 100. The blood circuit of thehemodialysis system 100 includes, among other things, a blood line setcomprising multiple blood lines 105, a drip chamber 106, and a dialyzer110. A blood pump (e.g., a peristaltic pump) 108 is configured to pumpblood through the blood circuit during treatment. The hemodialysissystem 100 also includes a dialysate circuit and various othercomponents that, for the sake of simplicity, are not described indetail. During hemodialysis treatment, blood is drawn from the patientand, after passing through the drip chamber 106, is pumped through thedialyzer 110 where toxins are removed from the blood and collected indialysate passing through the dialyzer. The cleansed blood is thenreturned to the patient, and the dialysate including the toxins(referred to as “spent dialysate”) is disposed of or recycled andreused. As discussed in greater detail below, during the hemodialysistreatment, drugs (such as erythropoiesis stimulating agents, e.g.,Epogen®, and iron replacement products, e.g., Venofer®) are alsodelivered to the drip chamber 106 using the drug delivery system 102.The drugs mix with the patient's blood within the drip chamber 106 andare then delivered to the patient along with the patient's blood.

As shown in FIG. 3, the modular drug delivery device 103 includes a drugvial holder 112 that defines four channels 114. Each of the channels 114is designed to hold captive a drug vial 116, 118. The channels 114 can,for example, be recesses that are formed within the drug vial holder 112and that are sized and shaped to receive only the caps and narrow neckportions of the vials 116, 118 such that the larger body portions of thevials sit above the holder 112. In the illustrated implementation, thevial 116 furthest to the left contains Venofer® and the three vials 118to the right of the Venofer® vial 116 contain Epogen®. Venofer® (ironsucrose injection, USP) is a sterile, aqueous complex of polynucleariron (III)-hydroxide in sucrose that is manufactured by American Regent,Inc. Venofer® is indicated in the treatment of iron deficiency anemia inpatients undergoing chronic hemodialysis who are receiving supplementalerythropoietin therapy. Epogen® is a drug that stimulates the productionof red blood cells and is also commonly used in dialysis patients.Epogen® is manufactured by Amgen, Inc.

The disposable drug administration fluid line set 109 is fluidlyconnected to each of the vials 116, 118. The drug administration fluidline set 109 includes four drug vial spikes 120 that connect to thevials 116, 118 in a manner to allow the drugs within the vials (e.g.,the Venofer® and Epogen®) to flow into feeder lines 122 via the drugvial spikes 120. Each of the feeder lines 122 is attached to aT-connector 124. The T-connectors 124 and associated tubing segments 126connect the feeder lines 124 to the drug delivery line 104. The drugvial spikes 120 can be formed of one or more relatively rigid medicalgrade plastics, such as polycarbonate or alphamethylstyrene (AMS), andthe various fluid lines can be formed of a more flexible medical gradeplastic, such as polyvinylchloride (PVC).

Each of the feeder lines 122, as shown in FIG. 3, passes through (e.g.,is threaded through) a bubble detector 128. The bubble detectors 128 arecapable of detecting air bubbles within the feeder lines 122. As aresult, each of the bubble detectors 128 can determine whether itsassociated drug vial 116, 118 is empty during treatment, because air isdrawn from the vial 116, 118 into the feeder line 122 when the vial isempty. In some implementations, the bubble detectors 122 are opticaldetectors. The OPB 350 bubble detector made by Optek can, for example,be used. Other types of optical detectors can alternatively oradditionally be used. Similarly, other types of sensors, such as sensorsutilizing ultrasound technology can be used as the bubble detectors.Examples of such sensors include the AD8/AD9 Integral UltrasonicAir-In-Line, Air Bubble Detector and the BD8/BD9 Integral Ultrasonic AirBubble, Air-In-Line & Liquid Level Detection Sensors (manufactured byIntrotek International (Edgewood, N.Y.)). In some implementations, thebubble detector 128 includes a sensor that, in addition to sensing thepresence of an air bubble within its associated feeder line 122, cansense the presence of the feeder line itself.

Downstream of the bubble detectors 128, the feeder lines 122 passthrough (e.g., are threaded through) occluders 130. Each of theoccluders 130 can be used to crimp the portion of the feeder line 122disposed therein to prevent fluid from passing through the feeder line122. In some implementations, the occluders 130 are solenoid based rams.Alternatively or additionally, other types of automated occluders can beused. The occluders 130 can be collectively operated in a manner suchthat only one feeder line 122 is unclamped at any particular time.

The drug delivery line 104 to which each of the feeder lines 122 isfluidly connected passes through (e.g., is threaded through) aperistaltic drug pump 132. The drug pump 132 includes multiple rollersthat compress the drug delivery line 104 in a manner to create a“pillow” of fluid (i.e., a “pillow” of air or liquid) that is pinchedbetween two points of the drug delivery line 104 that are compressed bythe pump rollers. The rollers are arranged around a circumference of arotatable frame. As the frame is rotated, the rollers force the “pillow”of fluid through the drug delivery line 104 toward the drip chamber 106(shown in FIG. 2). When the pump 132 is being operated and one of theoccluders 130 is open (i.e., not clamping its associated feeder line122), vacuum pressure is applied to the drug vial 116, 118 that isconnected to the feeder line 122 associated with the open occluder 130.In certain cases, the initial pressure in the drug vial 116, 118 isequal to the ambient pressure and when all of the drug has beendelivered, the ending pressure within the vial is about −10 psi. Inother words, the pressure within the drug vial 116, 118 progresses fromambient to −10 psi as the drug is delivered. The pump 132 is configuredto generate a vacuum pressure within the drug delivery line 104 andfeeder line 122 that exceeds the competing vacuum within the drug vial116, 118. As a result, the drug is drawn from the vial 116, 118, throughthe drug vial spike 120, through the feeder line 122, and into the drugdelivery line 104.

In some implementations, each channel 114 of the drug vial holder 112includes a sensor to sense the presence of a vial or drug container. Incertain implementations, each drug channel 114 includes a system whichidentifies the drug vial installed. The drug vial identification systemcan, for example, include a bar code reader that reads bar codes on thevials. Different types of sensors can alternatively or additionally beused. In some implementations, for example, the vial identificationsystem uses RFID technology. Other examples of suitable sensors includecolor sensors for sensing the color of color coded drug vials and/or forsensing the color of the drug within the vial, photo sensors (e.g.,cameras) that are equipped with text recognition software to read texton the drug vial, capacitive sensors that permit different size vials tobe detected, load cells or scales that detect the mass of the vial, andconductivity or electrical impedance sensors that can be used todetermine the type of drug within the vial.

The drug delivery device 103 also includes a control unit (e.g., amicroprocessor) that can power the various components of the drugdelivery device 103. The control unit can receive signals from and sendsignals to the various components of the drug delivery device 103,including, but not limited to, the bubble detectors 128, the occluders130, the drug pump 132, the drug vial ID sensors, and other sensorsalong the drug lines. The control unit can control the variouscomponents of the drug delivery device 103 based on information receivedfrom these components. For example, the control unit can control theoccluders 130 to ensure that only one of the occluders 130 is open at atime. This helps to ensure that drug is pulled from only one of thevials 116, 118 at a time during treatment. The control unit can alsodetermine the volume of drug delivered based on operation data of thedrug pump 132 and can control the occluders 130 based on the drug volumedetermined to have been delivered. For example, upon determining thatthe prescribed volume of the drug has been delivered, the control unitcan close the occluder 130 associated with that drug vial 116, 118 andopen the occluder 130 associated with the next drug to be delivered.

The control unit can also control the timing with which the variousoccluders 130 are opened and closed. For example, after the fullcontents of a vial have been evacuated, air will be sucked into thefeeder line 122 associated with that vial. As the air passes through thefeeder line 122, the bubble detector 128 will detect the air andtransmit a signal to the control unit indicating that the vial is empty.In response, the control unit can close the occluder 130 associated withthe empty vial and open the occluder 130 associated with the vialcontaining the next drug to be delivered. Upon receiving informationfrom the bubble detectors 128 indicating that all of the vials have beenemptied, the control unit can turn off the drug pump 132.

The control unit can also control certain components of the drugdelivery device 103 based on signals received from the drug vial IDsensors, which indicate the presence of a vial and/or the identity ofthe vial contents. Such an arrangement can help to ensure that thecorrect vials (e.g., the correct number of vials and the vialscontaining the correct contents) are used for the treatment. Uponreceiving signals from the drug vial ID sensors that do not match theinputted treatment information, for example, an alarm (e.g., an audibleand/or visual alarm) can be activated. Alternatively or additionally,the drug delivery device 103 can be configured so that treatment cannotbe initiated until the sensors detect the correct combination of vials.

The control unit can also be configured to adjust the patient'sundesired hematocrit and/or hemoglobin concentration to a value within adesired range at a predetermined time by i) employing a model forpredicting the patient's hematocrit and/or hemoglobin concentration at apredetermined time, the model including iron homeostasis, and by ii)operating the blood pump, the at least one ESA pump, and the at leastone iron replacement product pump to deliver blood, ESA, and ironreplacement product to the drip chamber via the at least one of theblood lines connected to the drip chamber and the at least two of thedrug lines connected to the drip chamber, according to an ESAadministration regimen predicted to adjust the patient's hematocritand/or hemoglobin concentration to a value within the desired range atthe predetermined time. Delivery of ESA and iron replacement product tothe drip chamber can be simultaneous, consecutive, or independent of oneanother.

The drug delivery device 103 (e.g., the control unit of the drugdelivery device 103) is configured to sense if the blood pump 108 of thedialysis machine 101 is running and to pause drug delivery if the bloodpump 108 is stopped. This technique prevents ‘pooling’ of the delivereddrug in the drip chamber 106 during treatment.

Still referring to FIG. 3, the drug delivery device 103 further includesa user interface 134 that is connected to the control unit. The userinterface 134 includes right/left arrow keys that allow the user tonavigate through displays associated with the vials 116, 118. The userinterface 134 also includes up/down arrow keys that enable the user toset the desired dosage for each of the vials 116, 118. In addition, theuser interface 134 includes start and stop keys that allow the user tostart and stop the drug delivery device 103.

Any of various other types of user interfaces can alternatively oradditionally be used. In some implementations, the drug delivery deviceincludes a user interface that allows the user to select a drug toinfuse from a menu. In certain implementations, the user may confirmthat the drug identified by the drug vial ID sensor is correct and/ormake appropriate adjustments. The user interface can be used to inputand/or monitor various different treatment parameters. Examples of suchparameters include drug dosage, drug delivery rate, amount of drugdelivered, status of the drug delivery for each drug channel, time,percent complete, percent remaining, time remaining, time delivered,date, patient ID, patient name, alarms, alerts, etc. Such userinterfaces can include a color graphical display. In certainimplementations, for example, the user interface is color codedaccording to drug, dosing, or status of drug delivery (e.g., done,running, ready, etc.).

The drug delivery device 103 also includes an alarm and/or alert systemto which the control unit of the drug delivery device 103 is connected.The alarm and/or alert system can be configured to emit a visual and/oraudio alarm and/or alert. The alarm and/or alert system can furtherinclude pre-programmed alarm and/or alert limitations so that when auser modifies any aspect of the system to be outside of the limitations,or the machine itself detects any aspects of the system to be outside ofthe limitations, the module emits an alarm and/or alert.

Turning back to modeling the life cycle of red blood cells, anerythrocyte (i.e., a red blood cell (RBC)) that reaches an age of 120days, has traveled approximately 480 kilometers, making about 170,000circuits, in each cycle enduring osmotic swelling and shrinkage and anumber of deformations while passing through capillaries. Theaccumulated damage to the membrane of the RBC is believed to lead to thedestruction of the cell.

About two decades ago it was observed that erythrocytes actually undergosuicidal death. The suicidal death of erythrocytes is called eryptosis.This term was chosen in order to point out the differences from andsimilarities to apoptosis, the programmed cell death of nucleated cells.Although there exist a number of studies concerning eryptosis, and someparts are understood in great detail, the precise process and mechanismsremain elusive.

Stimulation of eryptosis is governed by a complex signaling network andinvolves activation of cation channels. (F. Lang, C. Birka, S. Myssina,K. S. Lang, P. A. Lang, V. Tanneur, C. Duranton, T. Wieder, and S. M.Huber. Advances in Experimental Medicine and Biology, 559:211-217,2004). Shortly before they are phagocytosed by macrophages, which takesplace mainly in the spleen, one can observe distinctive cell shrinkage,plasma membrane microvesiculation and an exposure of phosphatidylserine(PS) on the cell surface. These PS-exposing erythrocytes are recognized,engulfed and degraded by macrophages. (D. Bratosin, J. Estaquier, F.Petit, D. Arnoult, B. Quatannens, J. P. Tissier, C. Slomianny, C.Sartiaux, C. Alonso, J. J. Huart, J. Montreuil, and J. C Ameisen. CellDeath Differ, 8(12):1143-1156, December 2001).

The normal lifespan and senescence can be either accelerated or delayedby environmental signals. Triggers of eryptosis include osmotic shock,oxidative stress, prostaglandin E2, energy depletion, chlorpromazine,aluminum, mercury, etc. Diseases which are associated with an earlyexpression of PS and thus accelerated eryptosis are, for instance, irondeficiency, sepsis, hemolytic uremic syndrome and sickle-cell anemia.See R. F. Zwaal, P. Comfurius, and E. B. Bevers Cellular and molecularlife sciences, 62:971-988, 2005. On the other hand eryptosis may beinhibited by erythropoietin, adenosine, catecholamines, nitric oxide andactivation of G-kinase. See M. Föller, S. M. Huber, and F. Lang. IUBMBLife, 60:661-668, 2008. Thus, EPO seems not only to increase the numberof circulating erythrocytes by preventing apoptosis of progenitor cellsbut also by prolonging the lifespan of circulating RBC by inhibiting thecation channels. See F. Lang, K. S. Lang, P. A. Lang, S. M. Huber, andT. Wieder Antioxidants & Redox Signaling, 8:1183-1192, 2006; A. B.Schwartz, B. Kelch, L. Terzian, J. Prior, K. E. Kim, E. Pequinot, and B.Kahn. ASAIO transactions, 36:M691M696, 1990. The cation channels arevolume-sensitive, and, after activation of the channel,phosphatidylserine asymmetry is reversed and phosphatidylserine isexposed to the outside of the erythrocyte. Binding of erythropoietin toRBC and inhibition of the channels may contribute to an increase in theerythrocyte number during erythropoietin therapy. See S. Myssina, S. M.Huber, C. Birka, P. A. Lang, K. S. Lang, B. Friedrich, T. Risler, T.Wieder, and F. Lang. Journal of American Society of Nephrology,14:2750-2757, 2003. Surprisingly, it was recently observed thaterythrocytes from erythropoietin overexpressing mice die faster ex vivo.See M. Foeller, R. S. Kasinathan, S. Koka, S. M. Huber, B. Schuler, J.Vogel, M. Gassmann, and F. Lang. American Journal ofPhysiology—Regulatory, Integrative and Comparative Physiology,293:R1127R1134, 2007. This observation leaves room for discussion. Ithas been suggested that: “Possibly, erythropoietin stimulates theexpression of genes in progenitor cells, which foster eryptosis and thusaccelerate erythrocyte death as soon as the erythropoietinconcentrations decline.” See F. Lang, E. Gulbins, H. Lerche, S. M.Huber, D. S. Kempe, and M. Föller. Eryptosis, Cellular Physiology andBiochemistry, 22:373-380, 2008.

Erythropoiesis and its Regulatory Mechanisms

The daily production rate in a healthy adult totals about 200×10⁹ redblood cells under normal conditions. This number can vary due to varyinginternal and environmental conditions. The body has to be able to adaptthe number of erythrocytes to situations of anemia and/or hypoxia (e.g.after bleeding, due to a shortened RBC lifespan, due to an enhancederyptosis, . . . ) as well as to situations of excessive red blood cells(e.g., high altitude dwellers descending to sea level).

The production of new erythrocytes takes place in the bone marrow.Undifferentiated stem cells in the bone marrow commit to the erythroidlineage and undergo a series of proliferations and differentiations. Theearliest stage of erythroid committed cells are called BFU-Es(Burst-Forming Unit Erythroid). Within a few days the BFU-Es mature toCFU-Es (Colony-Forming Unit Erythroid), then they undergo differentstages of erythroblasts, and finally they become reticulocytes. Thereticulocytes are released to the blood stream and within 1-2 days theyappear as mature erythrocytes. See Williams Hematology. Throughout thedescription below, the term progenitor cells is used to refer to BFU-Eand CFU-E cells, and the term precursor cells is used to subsumeproerythroblasts, basophilic erythroblasts, orthochromatophilicerythroblasts and bone marrow reticulocytes.

The primary control of erythropoiesis is governed by the hormoneerythropoietin. EPO is produced mainly by peritubular cells in the renalcortex based on a negative feedback mechanism. The kidneys detect thepartial pressure of oxygen in blood and react by releasing more EPO ifthe oxygen content is too low and vice versa. The concentration of EPOaffects the number of circulating RBC by determining the number of cellsthat mature into erythrocytes. Erythropoietin plays a role in therecruitment of stem cells to the erythroid lineage, prevents apoptosisof progenitor cells and affects the maturation velocity of progenitorand precursor cells. Thus, disturbances in oxygen delivery can beadjusted for by an adaptive resetting of the rate of erythropoiesis.Additionally, there exists a physiological process wherein macrophagesstart to phagocytose young erythrocytes (neocytes), which is calledneocytolysis. Neocytolysis seems to be triggered by a drop in the EPOlevel and helps to regulate situations of excessive red cell mass. SeeRice 2001. Further, recent studies suggest that the concentration of EPOin blood influences the clearance of senescent RBCs and that EPO canprolong the lifespan of erythrocytes by inhibiting cation channels.Progenitor and Precursor Cells

When a stem cell is committed to the erythroid lineage, it develops intoan erythrocyte. The number of stem cells which enter the differenthematopoietic lineages is determined by interleukines and growthfactors. The immature erythrocyte undergoes a number of changes instructure, appearance and its requirements during the process ofdifferentiation and proliferation. BFU-E cells, the first committederythroid cells, express only a very small number of EPO-receptors(EpoR) and therefore they are almost EPO-independent for their survival.Within a few days they develop into CFU-E cells.

As cells mature to CFU-Es, the number of EpoR on the cell surfaceincreases distinctly, and the cells become absolutely dependent onerythropoietin to prevent them from apoptosis. See H. Wu, X. Liu, R.Jaenisch, and H. F. Lodish. Cell, 83(1):59-67, 1995 (hereinafter “Wu etal.”). Under normal conditions, large numbers of generated CFU-E are notsurviving. After that very EPO-sensitive phase, the density of EpoRsdeclines sharply on early erythroblasts and EpoRs almost disappear atthe stage of orthochromatophilic erythroblasts. See J. P. Greer, J.Foerster, G. M. Rodgers, F. Paraskevas, B. Glader, D. A. Arber, and R.T. Jr. Means. Wintrobe's Clinical Hematology, volume 1. LippincottWilliams & Wilkins, 12th edition, 2009 (hereinafter “Wintrobe's ClinicalHematology”). Under high levels of EPO, the marrow transit time ofprecursor cells is shortened, and high EPO levels result in release ofstill immature reticulocytes, which are referred to as stressreticulocytes. (See Williams Hematology).

One significant difference between progenitor and precursor cells is thesynthesis of hemoglobin. Unlike progenitor cells, all types of precursorcells synthesize heme and the increased demand for iron for this processis reflected by a sharp increase in the number of transferrin receptors(TfR) in proerythroblasts (transferrin receptors are the path of ironuptake). TfRs reach their peak in intermediate erythroblasts and declineafterward and mature bone marrow reticulocytes express only a few TfRs(in comparison). A diagram of the different erythroid cell types isshown in FIG. 1.

Neocytolysis

The body has to be able to deal with situations of too few circulatingerythrocytes as well as with excessive red blood cells. Whereas it isvery effective to change the rate of erythropoiesis to adapt the systemto situations when red cell numbers are low, this is not an adequateregulatory mechanism in situations of red cell mass excess. This isbecause a decrease in the rate of erythropoiesis takes relatively longto have a noticeable effect on the number of circulating red bloodcells, due to the long lifespan of erythrocytes.

A mere 20 years ago, preventing apoptosis of red cell progenitors hadbeen thought to be the sole regulating effect of EPO. Because this wouldonly allow for a slow adaption in situations of too many RBC, thecontrol of the body of red cell mass would be coarse. A decline inerythropoietin level results in more progenitor cells dying, but asuppression of the hormone has no effect on maturation of erythroidprecursors. Thus, no decrement in the erythrocyte production would beobservable within one week after EPO has declined. However, studies doneon residents at very high altitude who rapidly descended to sea level,showed a decrease in red cell mass of 10-18% in the first 7-10 days.(See Rice 2001). Persons living in hypoxic environments, like highaltitude, are normally polycythemic. Thus, their hematocrit is too highunder normal conditions and as a consequence the release of EPO issuppressed.

Further investigations of situations where EPO levels are lower thannormal (polycythemic high altitude dwellers, anemia of renal failure,human model based on EPO administration then withdrawal), suggest that asuppression of EPO leads to selective hemolysis of young red bloodcells. (See Rice 2005). This process is called neocytolysis, to stressthe fact that young erythrocytes (i.e., neocytes) are uniquelysusceptible. See C.-C. Chang, Y. Chen, K. Modi, O. G. Awar, C. P.Alfrey, and L. Rice. Journal of Investigative Medicine, 57:650-654,2009; Rice 2005; M. M. Udden, Pickett M. H. Driscoll, T. B., C. S.Leach-Huntoon, and C. P. Alfrey. The Journal of Laboratory and ClinicalMedicine, 125:442-449, 1995.

Neocytolysis is initiated by a fall in erythropoietin levels and in Rice2005 it was shown that, for instance, low doses of EPO, administered tohigh altitude dwellers on descent, completely abrogated neocytolysis. Atthe moment it remains unclear whether it is the rate of decline in EPOlevel or the drop of EPO beneath a certain threshold that acts as atrigger for neocytolysis. See C. P. Alfrey and S. Fishbane. NephronClinical Practice, 106:149-156, 2007 (hereinafter “Alfrey et al.”); Rice1999. The current approach to treatment of anemia with ESAs deviatessignificantly from the normal internal systemic process. Neocytolysis,for instance, contributes to renal anemia. It is precipitated by thepathologic endogenous erythropoietin deficiency of renal disease and theshort bursts in EPO concentration followed by sharp declines in serumlevels due to administration of the hormone (especially duringintravenous (i.v.) administration). It would be desirable to choosedosing schedules such that neocytolysis is minimized or totallyprevented. (See Alfrey et al.; Rice 1999). A mathematical model could bevery useful to define better ESA administration schemes which stimulateerythropoiesis in a more natural way and abrogate neocytolysis. Betterand/or alternative administration schemes can include daily ESAadministration, continuous ESA administration (e.g., via a pump), ormonthly ESA administration.

Iron and Erythropoiesis

Erythropoiesis is a very complex process and even though erythropoietinis the key regulator, there are other proteins (e.g., interleukins,etc.) and substances (e.g., iron, folic acid, vitamin B₁₂, etc.) thatare needed for an optimal erythropoiesis. Iron is a very critical factorfor an effective production of red blood cells. The metal isindispensable for hemoglobin synthesis and the hemoglobin protein is theactual oxygen transporter in erythrocytes. The molecule makes up about97% of the dry weight of red blood cells. If the body is not able toprovide sufficient iron for the production of erythrocytes, thenimpaired erythropoiesis will result. (S. Finch, D. Haskins, and C. A.Finch. The Journal of Clinical Investigation, 29:1078-1086, 1950;Williams Hematology). In general, a healthy adult has difficultyproviding sufficient iron to support production rates greater than threetimes basal. Higher rates may be possible when administering EPO andiron and in some diseases, e.g. hemochromatosis. (L. T. Goodnough.Nephrology Dialysis Transplantation, 17:14-18, 2002 (hereinafter“Goodnough 2002”); C. A. Finch. Blood. The Journal of American Societyof Hematology, 60(6):1241-1246, 1982) (hereinafter “Finch 1982”). Anundersupply of iron additionally leads to an increase in the number ofhypochromic RBC. Hypochromic cells (i.e., cells in which the amount ofhemoglobin is lower than the normal 26 pg) are small, have a reducedoxygen carrying capacity and are relatively fragile, and thus are likelyto die earlier than normochromic erythrocytes. See D. M. Wrighting andN. C. Andrews. Iron Homeostasis and erythropoiesis. Current Topics inDevelopmental Biology, 82:141-159, 2008 (hereinafter “Wrighting”); A.Loria, L. Sánchez-Medal, R. Lisker, E. De Rodríguez, and J. Labardini.Br J Haematol, 13(3):294-302, May 1967.

Erythroid precursor cells are the most avid consumers of iron in thebody. See P. Ponka and C. N. Lok. The International Journal ofBiochemistry & Cell Biology, 31:1111-1137, 1999 (hereinafter “Ponka1999”). The immature erythrocyte has only a few days to synthesize allthe hemoglobin which the mature cell contains. Hemoglobin is ametalloprotein. The name refers to its special structure. A hemoglobinmolecule consists of four subunits each composed of a globular proteinembedding a heme group, and every heme group in turn contains one ironatom. Erythroid cells rely completely on transferrin receptors to takeup iron. See R. Y. Chan, C. Seiser, H. M. Schulman, L. C. Kuhn, and P.Ponka. European Journal of Biochemistry, 220:683-692, 1994. Unlikeprogenitor cells, all types of precursor cells synthesize heme and theincreased demand for iron for this process is reflected by a sharpincrease in the number of TfR in proerythroblasts. Transferrin receptorsreach their peak in intermediate erythroblasts followed by a decreasewith further maturation.

In precursor cells, heme is essential for maintaining a “normal” numberof TfR. Studies showed that inhibition of heme synthesis stronglyinhibited TfR expression. It seems that there exists a positive feedbackmechanism in which heme promotes a high rate of transferrin receptorsynthesis. A high number of transferrin receptors enhances iron uptakeand that in turn keeps hemoglobin synthesis at high levels. See P.Ponka. Blood, 89:1-25, 1997. There is evidence that precursor cellswhich have a low hemoglobin content are still able to proliferate but donot differentiate and undergo apoptosis. See J. A. Schmidt, J. Marshall,M. J. Hayman, P. Ponka, and H. Beug. Cell, 46:41-51, 1986 (hereinafter“Schmidt”).

1. Iron Homeostasis

Iron plays a vital role in metabolic processes in all living organisms,including oxygen transport, DNA synthesis, and electron transport. It isthe most important essential trace element and amounts to about 35-45and 50-55 mg/kg body weight in adult women and men, respectively. SeeWilliams Hematology; Ponka 1999. For a very good overview of ironhomeostasis and metabolism, see Crichton and Lieu. Homeostasis is thestate of equilibrium (balance between opposing pressures) in the bodywith respect to various functions and to the chemical compositions ofthe fluids and tissues. See Steadman's Medical Dictionary, 26^(th)Edition, Williams & Wilkins, 1995.

Iron is a key component of many cellular enzymes (e.g. oxidases,catalases, peroxidases, cytochromes, etc.). These enzymes are criticalin many basic cellular processes, including cell proliferation anddifferentiation, DNA and RNA synthesis and electron transport. Further,iron, bound in heme, is essential for oxygen transport. Under normalconditions, total body iron burden is very tightly regulated, becauseexcessive iron leads to tissue damage. Once an atom of iron enters thebody, it is conserved with remarkable efficiency and can remain in thebody for more than ten years. Iron is only lost via loss of cells(epithelial cells, bleeding) and very small amounts are lost via urineand sweat. Unlike other mammals, humans are not able to actively excreteiron. Thus, because of the toxic nature of iron, it is very important tokeep absorption and needs well balanced.

Iron is absorbed from dietary sources via the duodenum. The uptake isregulated by the hormone hepcidin. Under normal conditions, the dailyabsorption of iron amounts to 1-2 mg. In general, in the organism ironis bound to proteins, because free iron results in formation of freeradicals and destroys cells. The iron exporter ferroportin is needed toegress the metal from cells (e.g. enterocytes, macrophages,hepatocytes). Ferroportin is not only the effector of cellular ironexport but also the receptor of the iron regulating hormone hepcidin.See Williams Hematology.

Much of the iron in the human body can be found in circulatingerythrocytes incorporated in hemoglobin (1 ml of packed cells containabout 1 mg of iron). However, most of the iron used for the productionof new erythrocytes comes from hemoglobin recycling and not fromabsorption. When a senescent red cell is phagocytosed, the macrophageinternalizes hemoglobin, degrades it, and exports it back out into bloodcirculation.

Iron is stored in the form of ferritin or hemosiderin, where the latteris more a kind of a long-term storage, because iron within deposits ofhemosiderin is very poorly available. Storage sites can be found mainlyin the liver, the spleen, and the bone marrow but, for instance, smallamounts of hemosiderin can be found in almost any tissue. There aredifferences between men and women in the amount of iron that is stored(men: about 200-400 mg, women: about 1000 mg,). See C. Crepaldi, A.Brendolan, V. Bordoni, M. R. Carta, V. D. Intini, F. Gastaldon, P.Inguaggiato, and C. Ronco. Hemodialysis International, 7(3):216-221,2003 (hereinafter “Crepaldi”).

Iron is carried around in plasma by a transport protein, which is calledtransferrin. Transferrin is produced by the liver. Cells which are inneed for iron express transferrin receptors, and, after transferrinforms a complex with this receptor, iron is transported into the cell.Transferrin exists in three different forms: iron-free (apotransferrin),one iron (monoferric transferrin) and two iron (differic transferrin)atoms bound. The concentration of iron and the amount of transferrinpresent in plasma determine the relative abundance of each form.Although iron bound to transferrin amounts to only about 3 mg of totalbody iron, it is the most important iron pool with the highest rate ofturnover, about 30 mg/day. Approximately 80% of this iron is transportedto the bone marrow and used by erythroid cells. See Lieu. The amount oftransferrin-bound iron is determined by three processes: macrophage ironrecycling, iron storage (hepatic) and duodenal iron absorption. SeeWrighting.

In the last 15 years, understanding of the regulation of systemic ironhomeostasis has changed substantially. The antimicrobial peptidehepcidin (hepatic bactericidal protein), which is secreted by the liver,was identified to be the systemic iron regulator.

Hepcidin inhibits iron transport by binding to the iron exporterferroportin. Hence, it keeps gut enterocytes from secreting iron intocirculation, thereby functionally reducing iron absorption. Further, itprevents iron release from macrophages and blocks the stores as well, byshutting off the means of transport out of cells. See FIG. 4 for adiagram of the iron cycle and daily turnover rates and FIG. 5 for asketch of systemic iron regulation.

There exist some sensing mechanisms in the liver that gauge the amountof iron needed to support erythropoiesis. Note that hepatocytes in theliver, which secrete hepcidin, also release transferrin into circulationand are one of the major iron storage sites. If the concentration ofiron in stores is too low and/or erythropoiesis is increased, then thelevel of the hormone hepcidin is decreased. Thus, for instance, lowhepcidin concentrations are observed in patients with absoluteiron-deficiency anemia, or anemias with high erythropoietic activity.See L. T. Goodnough, E. Nemeth, and T. Ganz. Blood, 116:4754-4761, 2010(hereinafter “Goodnough 2010”). The consequence of a lower hepcidinlevel is that more iron is taken up via the duodenum and more iron isreleased from macrophages and from the stores. See Crichton; Fleming. Ifiron stores are full and/or erythropoiesis is decreased, then morehepcidin is secreted by hepatocytes. As a result of the binding to theiron exporter ferroportin, dietary iron absorption, macrophage ironrecycling, and release of iron from stores are partly or fully blocked.

Hepcidin is an acute phase-reactant, and, during inflammation, cytokinesstimulate overproduction of hepcidin. Hence, the feedback is alteredduring inflammation, which is a prominent problem in chronic kidneydisease. As a consequence, the body is not able to sufficiently supplydeveloping red blood cells with iron and erythropoiesis is impaired,even though there might be more than enough iron in stores andmacrophages, but it is simply not available. This paradoxical situationis called functional iron deficiency, whereas a situation when the ironcontent of the body is too low (depleted stores) is referred to asabsolute iron deficiency.

2. A Mathematical Model for Erythropoiesis Including Iron Homeostasis

For healthy persons, one can assume that erythropoiesis is alwayssufficiently supplied with iron. This assumption certainly does not holdfor some types of anemia and it is not valid for the majority ofdialysis patients. Therefore, the model developed in this sectionincorporates an extensive iron model in order to expand the possibleapplications to all dialysis patients and to pathologies where anunbalanced iron homeostasis leads to an abnormal erythropoiesis (e.g.,anemia of absolute and functional iron deficiency, hemochromatosis,etc.). The incorporation of an iron submodel and its effects onerythropoiesis has a significant influence on the structure of thepopulation equations.

A reduced reaction of the body to erythropoietin or ESAs is often theresult of an insufficient iron availability in the body. A model whichgathers the balance between developing precursor cells and iron that canbe delivered to the bone marrow for hemoglobin synthesis is needed.Further, iron homeostasis itself is a very complex process. There arecomplex feedback mechanisms which determine how much iron is shiftedaround from one place to another in the body and how much is absorbedvia dietary input. Unfortunately, most of these feedback mechanisms arenot well understood.

A compartmental model for iron is considered and linked to a structuredpopulation model which governs the dynamics of the erythroid cellpopulation. For a detailed explanation of amendments and modificationsin the population equations for the iron model see Section 2.1. InSection 2.2 the derivation of the iron model is explained. For anoverview of how the model developed in this section is organized, seeFIG. 6B.

2.1. Population Model.

The structured population equations are designed to link the iron modelto the cell population model in a physiologically meaningful way. Fivedifferent cell population classes (BFU-E, CFU-E, erythroblasts, marrowreticulocytes, circulating RBCs) are considered. The inclusion of ironmakes it necessary to consider some physiological processes in detail.In many applications of structured population equations, the onlyattribute used to distinguish cells is the cell age. This concept isinsufficient for the population class of the precursor cells, where ironis used to synthesize hemoglobin. Here, cells are characterized by theirhemoglobin content, and the number of tranferrin receptors on the cellmembrane, in addition to the cell age. For mature erythrocytes, cell ageand hemoglobin content are used as attributes.

The two equations describing the progenitor cell populations (BFU-Es andCFU-Es) are solely structured using cell age. These types of cells donot synthesize hemoglobin and, therefore, they are not among the majorconsumers of iron. Hence, these two classes can be described bypopulation equations of the form

${{{\frac{\partial}{\partial t}{u\left( {t,\mu} \right)}} + {{v\left( {E(t)} \right)}\frac{\partial}{\partial\mu}{u\left( {t,\mu} \right)}}} = {\left( {{\beta\left( {{E(t)},\mu} \right)} - {\alpha\left( {{E(t)},\mu} \right)}} \right){u\left( {t,\mu} \right)}}},$where u(t,μ) is the population density of cells at time t with cell ageμ. The functions β and α describe the rate of division and rate ofapoptosis, respectively, of the cells. Further, v describes thematuration velocity.

Precursor cells are the only erythroid cells that synthesize hemoglobin.In order to introduce these phenomena in the model and to keep track ofthe amount of hemoglobin the cells are synthesizing, two additionalattributes are introduced in addition to cell age:

-   -   the number τ of TfRs per cell and    -   the amount h of hemoglobin each cell is carrying.        Thus, the general population equation for precursor cells is

$\begin{matrix}{{{\frac{\partial}{\partial t}{u\left( {t,\mu,h,\tau} \right)}} + {\frac{\partial}{\partial\mu}\left( {{v_{1}\left( {t,\mu,h,\tau} \right)}{u\left( {t,\mu,h,\tau} \right)}} \right)} + {\frac{\partial}{\partial h}\left( {{v_{2}\left( {t,\mu,h,\tau} \right)}{u\left( {t,\mu,h,\tau} \right)}} \right)} + {\frac{\partial}{\partial\tau}\left( {{v_{3}\left( {t,\mu,h,\tau} \right)}{u\left( {t,\mu,h,\tau} \right)}} \right)}} = {{8\;{\beta\left( {t,\mu,{2h},{2\tau}} \right)}{u\left( {t,\mu,{2h},{2\tau}} \right)}} - {\left( {{\beta\left( {t,\mu,h,\tau} \right)} + {\alpha\left( {t,\mu,h,\tau} \right)}} \right){{u\left( {t,\mu,h,\tau} \right)}.}}}} & (1)\end{matrix}$Here u(t,μ,h,τ) denotes the population density of the cell population attime t with maturity μ, hemoglobin content h and number of transferrinreceptors τ. Further, the functions β and α describe the rate ofdivision and rate of apoptosis of the cells. The functions v₁, v₂, v₃define the rate with which the structural variables μ, h and τ arechanging, respectively. For a derivation of the partial differentialequation (1), see below.

Circulating red blood cells carry hemoglobin but they do not expressTfRs for synthesis of the hemoglobin molecule and, therefore, adistinction is made between cells using only cell age μ and hemoglobincontent h as structural variables. Hence, the general populationequation reads as follows

${{\frac{\partial}{\partial t}{u\left( {t,\mu,h} \right)}} + {\frac{\partial}{\partial\mu}\left( {{v_{1}\left( {t,\mu,h} \right)}{u\left( {t,\mu,h} \right)}} \right)} + {\frac{\partial}{\partial h}\left( {{v_{2}\left( {t,\mu,h} \right)}{u\left( {t,\mu,h} \right)}} \right)}} = {{- {\alpha\left( {t,\mu,h} \right)}}{{u\left( {t,\mu,h} \right)}.}}$2.1.1. Progenitor Cells: BFU-E and CFU-E Cells.Assumptions:

-   -   1. The number of stem cells, which commit to the erythroid        lineage depends on EPO.    -   2. Cells normally stay in this stage for 8 days (3 days BFU-E        and 5 days CFU-E).    -   3. The number of transferrin receptors on BFU-E and CFU-E cells        is negligible.    -   4. EPO has no effect on the number of divisions or the rate of        apoptosis of BFU-E.    -   5. The proliferation rate of CFU-E cells is constant.    -   6. The rate of apoptosis depends strongly on EPO levels.    -   7. The maturation velocities of BFU-E and CFU-E cells are        constant.

Different kinds of growth factors affect the number of stem cells and invitro studies suggest that high levels of EPO can result in up to 20%more stem cells committing to the erythroid lineage. BFU-E and CFU-Ecells express only a very small number of transferrin receptors on thecell membrane and the influence of iron on these cells is neglected.This is not entirely true, because, like all proliferating cells, theyneed a small amount of iron for division. These amounts are so smallcompared to the consumption of iron for hemoglobin synthesis that theyare neglected. Even in situations of iron deficiency, the assumption isthat there is enough iron for cell divisions of progenitor cellsavailable and only hemoglobin synthesis is impaired. Thus, the followingpopulation equation is obtained for BFU-E cells:

$\quad\left\{ \begin{matrix}{{{{\frac{\partial}{\partial t}{p\left( {t,\mu^{p}} \right)}} + {\frac{\partial}{\partial\mu^{p}}{p\left( {t,\mu^{p}} \right)}}} = {\beta^{p}{p\left( {t,\mu^{p}} \right)}}},} \\{{{p\left( {t,0} \right)} = {S_{0}\left( {E(t)} \right)}},} \\{{{p\left( {0,\mu^{p}} \right)} = {p_{0}\left( \mu^{p} \right)}},}\end{matrix} \right.$where p(t,μ^(p)) is the population density of the cell class at time twith maturity μ^(p), 0≦μ^(p)≦μ_(max) ^(p)=3, t>0. Further, β^(p) is theconstant proliferation rate (Assumption 4) and p₀(μ^(p)) is thepopulation density at t=0. The number S₀(E(t)) of cells committing tothe erythroid lineage is given by a sigmoidal function,

${{S_{0}\left( {E(t)} \right)} = {\frac{a_{1} - b_{1}}{1 + {\mathbb{e}}^{{{- k_{1}}{E{(t)}}} + c_{1}}} + b_{1}}},$where E(t) is the EPO concentration at time t and a₁, b₁, c₁ and k₁ arepositive constants and a₁>b₁. The function S₀ increases monotonicallywith increasing EPO concentration (Assumption 1).

The equation for CFU-E cells is:

$\quad\left\{ \begin{matrix}{{{{\frac{\partial}{\partial t}{q\left( {t,\mu^{q}} \right)}} + {\frac{\partial}{\partial\mu^{q}}{q\left( {t,\mu^{q}} \right)}}} = {\left( {\beta^{q} - {\alpha^{q}\left( {E(t)} \right)}} \right){q\left( {t,\mu^{q}} \right)}}},} \\{{{q\left( {t,\mu_{m\; i\; n}^{q}} \right)} = {p\left( {t,\mu_{{ma}\; x}^{p}} \right)}},} \\{{{q\left( {0,\mu^{q}} \right)} = {q_{0}\left( \mu^{q} \right)}},}\end{matrix} \right.$where q(t,μ^(q)) is the population density of the CFU-E class at time twith maturity μ^(q), t>0 and 3=μ_(min) ^(q)≦μ^(q)≦μ_(max) ^(q)=8. Theconstant β^(q) describes the division rate (Assumption 4) andα^(q)(E(t)) the apoptosis rate which depends on the EPO-concentration(Assumption 5). Moreover, the number of cells leaving the BFU-E cellstage and entering the CFU-E cell stage per unit time requires theboundary condition q(t,μ_(min) ^(q))=p(t,μ_(max) ^(p)). Finally,q₀(μ^(q)) is the population density at t=0.

Again, a sigmoidal function is used to describe the dependence of therate of apoptosis on EPO,

${{\alpha^{q}\left( {E(t)} \right)} = {\frac{a_{2} - b_{2}}{1 + {\mathbb{e}}^{{k_{2}{E{(t)}}} - c_{2}}} + b_{2}}},$where E(t) is the EPO concentration and a₂, b₂, c₂ and k₂ are positiveconstants, a₂>b₂. The function α^(q) is per definition a monotonicallydecreasing function. Thus, a higher level of EPO causes more cells tosurvive.

Finally, because of Assumption 7, it is assumed that cell age equals theactual age of the cells, i.e., v^(p)=v^(q)≡1, ∀t≧0.

2.1.2. Precursor Cells: Erythroblasts and Marrow Reticulocytes.

Assumptions:

-   -   8. The class erythroblasts consists of all cell stages from        proerythroblast to orthochromatophilic erythroblast.    -   9. Cells normally stay in this stage for 6-8 days (5 days        erythroblasts and 1-3 days marrow reticulocytes).    -   10. Under high levels of EPO, so called “stress reticulocytes”        are released.    -   11. All types of precursor cells synthesize hemoglobin.    -   12. The rate at which the amount of hemoglobin in precursor        cells changes, depends on the number of TfRs on the cell surface        and the currently available amount of iron.    -   13. The rate of change of TfRs in precursor cells depends on the        cell age and the hemoglobin content of the cell.    -   14. Cells entering the erythroblast class may be equipped with        different numbers of TfRs, i.e., TfRs need not be uniformly        distributed among the erythroblast population at time t=0.    -   15. Hemoglobin is conserved during cell division and a mother        cell divides into 2 daughter cells carrying equal amounts of        hemoglobin.    -   16. TfRs are conserved during division and a mother cell divides        into 2 daughter cells with equal numbers of TfRs.    -   17. The maturation velocity of erythroblasts is constant.    -   18. The proliferation rate of erythroblasts is constant.    -   19. The rate of apoptosis of erythroblasts depends on cell age        and hemoglobin.    -   20. The maturation velocity of reticulocytes depends on EPO.    -   21. The rate of apoptosis of reticulocytes depends on cell age        and on hemoglobin.    -   22. Reticulocytes mature but do not proliferate.    -   23. Iron that enters a precursor cell is incorporated into        hemoglobin, i.e., no free iron in the cell is considered.

The number of TfRs increases sharply in the proerythroblasts and reachesits peak in intermediate erythroblasts before declining again(proerythroblasts: 300,000 TfRs/cell; basophilic erythroblasts: 800,000TfRs/cell; mature reticulocytes: 100,000 TfRs/cell). See Wintrobe'sClinical Hematology. This reflects the increased demand of the cells foriron to synthesize hemoglobin Almost all iron contained in the cells issynthesized to hemoglobin, leaving mature erythrocytes with only anegligible amount of non-heme iron (Assumption 23). Depending on thematurity of cells different numbers of TfRs are to be expected(Assumption 13). Cells that are saturated with hemoglobin reduce thelevel of receptors. Further, there is evidence that in contrast to othercells, in the precursor cells the number of transferrin receptors isalso controlled by a positive feedback involving the amount of alreadysynthesized hemoglobin in the cell. In vitro hemoglobin was shown to beessential for these cells to maintain a normal number of TfRs. See Ponka1999.

For the erythroblasts population class, the following PDE is obtained

$\quad\left\{ \begin{matrix}\begin{matrix}{{\frac{\partial}{\partial t}{r\left( {t,\mu^{r},h^{r},\tau^{r}} \right)}} + {\frac{\partial}{\partial\mu^{r}}r\left( {t,\mu^{r},h^{r},\tau^{r}} \right)} +} \\{{{v_{2}^{r}\left( {\tau,{a_{1}(t)}} \right)}\frac{\partial}{\partial h^{r}}{r\left( {t,\mu^{r},h^{r},\tau^{r}} \right)}},}\end{matrix} \\\begin{matrix}{{{+ \frac{\partial}{\partial\tau^{r}}}\left( {{v_{3}^{r}\left( {\mu^{r},h^{r},\tau^{r}} \right)}{r\left( {t,\mu^{r},h^{r},\tau^{r}} \right)}} \right)} = {{8\;\beta\; u\left( {t,\mu^{r},{2h^{r}},{2\tau^{r}}} \right)} -}} \\{{\left( {\beta + {\alpha^{r}\left( {\mu^{r},h^{r}} \right)}} \right){r\left( {t,\mu^{r},h^{r},\tau^{r}} \right)}},}\end{matrix} \\{{{r\left( {t,\mu_{m\; i\; n}^{r},0,\tau^{r}} \right)} = {{q\left( {t,\mu_{{ma}\; x}} \right)} \cdot {\eta\left( \tau^{r} \right)}}},} \\\begin{matrix}{{{r\left( {t,\mu^{r},h^{r},\tau^{r}} \right)} = 0},} & {{if}\mspace{14mu}\left( {\mu^{r},h^{r},\tau^{r}} \right)\mspace{14mu}{is}\mspace{14mu} a\mspace{14mu}{boundary}\mspace{14mu}{point}} \\\; & {{{{and}\mspace{14mu} h^{r}} > 0},}\end{matrix} \\{{{r\left( {0,\mu^{r},h^{r},\tau^{r}} \right)} = {r_{0}\left( {\mu^{r},h^{r},\tau^{r}} \right)}},}\end{matrix} \right.$where r(t,μ^(r),h^(r),τ^(r)) is the population density of theerythroblast class at time t with maturity μ^(r), hemoglobin contenth^(r) and number of TfRs τ^(r),8=μ_(min) ^(r)≦μ^(r)≦μ_(max) ^(r)=13, 0≦h^(r)≦h_(max) ^(r), τ_(min)^(r)≦τ^(r)≦τ_(max) ^(r), t>0. The population density at t=0 is given asr₀(μ^(r),h^(r),τ^(r)). Assumption 17 allows setting v_(μ) ^(r)≡1, i.e.,the cell age of the cell coincides with the actual age of the cell. Thecontact rate between the TfRs of a cell and transferrin molecules inplasma carrying iron is governed by the law of mass action:d ₁τ^(r)(t)α₁(t),where d₁ is a rate constant, τ^(r)(t) is the number of TfRs on the cellmembrane at time t and a₁(t) is the concentration of iron in plasma attime t. Thus, the amount of hemoglobin in a cell evolves according tov ₂ ^(r)(τ^(r) ,a ₁)=d ₂τ^(r) a ₁,  (2)where d₂ is a rate constant that is different from d₁ (note that fouratoms of iron are needed for one hemoglobin molecule). The change of thenumber of TfRs on the cell surface is described by a function v_(τ)^(r), depending on the cell age μ^(r) and the amount of hemoglobin h^(r)(Assumption 13). The function v_(τ) ^(r) is chosen such that itincreases in the beginning and then decreases when the cells are “fullysaturated” with hemoglobin.v ₃ ^(r)(μ^(r) ,h ^(r),τ^(r))=a ₃τ^(r)(h _(*) ^(r) −b ₃(h ^(r))²)−c₃max(τ^(r),0)(μ^(r)−μ_(max) ^(q))/(1+h ^(r)),where a₃, h_(*) ^(r), b₃ and c₃ are constants and μ_(max) ^(q) is themaximal maturity of progenitor cells. The second term assures that cellswhich progress in maturation, but do not accumulate sufficient iron,lose transferrin receptors. The term max(τ^(r),0) is needed to guaranteethat the number of TfRs can not become negative. The maximum function isnot differentiable, therefore it is approximated by the followingfunctionsmax(τ^(r))=0.5τ^(r)+0.5√{square root over ((τ^(r))²+0.01)}.

The rate of division β^(r) is constant (Assumption 15). Note that theterm 8βu(t,μ^(r),2h^(r),2τ^(r)) in the equation arises because ofAssumptions 15 and 16. For further details see below. To model the rateof apoptosis α^(r), a decrease in expression of TfRs is taken intoaccount and leads to a decrease in iron uptake and, thus, lesshemoglobin is synthesized. Cells containing small amounts of hemoglobinarrest in differentiation (proliferation is not affected) and this leadsto premature cell death. See Schmidt Therefore, it is assumed that α^(r)depends on μ^(r) and h^(r) (Assumption 19)α^(r)(μ^(r) ,h ^(r))=b ₄(h**−h ^(r))²/(1+(μ**−(μ^(r)−μ_(max)^(q)))²),  (3)where b₄, h**, μ**>0 are constants. The parameters h** and μ** can, forinstance, be chosen ash**=h _(max) ^(r), μ**=μ_(max) ^(r)−μ_(max) ^(q),  (4)where h_(max) ^(r) is the maximal hemoglobin amount of precursor cellsand μ_(max) ^(r) is the maximal maturity of precursor cells.

In principle, the boundary conditions arise as a result of the fact thatcells from the precedent class leave the class when they reach themaximum maturity and enter the next population class. Hence, there is aconstant flux of cells from the CFU-E class to the erythroblasts class.However, the passage is more complicated because there is a flux from apopulation class with only one structural variable to a population classwith three attributes. Thus, it has to carefully considered at whichpoints of the boundary the cells are entering. It is useful to keep thephysiological situation in mind. Progenitor cells do not synthesizehemoglobin. Hence, there are no cells entering the erythroblastspopulation class with h>0. On the other hand, it is not natural to applythe same assumption for transferrin receptors. For although the numberof TfRs on progenitor cells is neglected, in fact, very mature CFU-Eexpress a number of TfRs. Further, it seems reasonable to assume that atthe time when they become proerythroblasts not all cells incorporateexactly the same amount of TfRs on the surface, but that there is somedistribution of this attribute among the cells (Assumption 14). Thus,the concept for the boundary includes, that∫_(τ) _(min) ^(τ) ^(max) η(τ)dτ=1,where, for instance, η(τ) can be chosen such that the structuralvariable τ is normally distributed with a certain mean value andvariance. For an illustration of this situation see FIG. 7.

The revisited population equation for marrow reticulocytes reads

$\quad\left\{ \begin{matrix}\begin{matrix}{\left. {\frac{\partial}{\partial t}{s\left( {t,\mu^{s},h^{s},\tau^{s}} \right)}} \right) + {{v_{1}^{s}\left( {E(t)} \right)}\frac{\partial}{\partial\mu^{s}}{s\left( {t,\mu^{s},h^{s},\tau^{s}} \right)}} +} \\{\left. {{v_{2}^{s}\left( {\tau^{s},{a_{1}(t)}} \right)}\frac{\partial}{\partial h^{s}}{s\left( {t,\mu^{s},h^{s},\tau^{s}} \right)}} \right),}\end{matrix} \\{{{{+ \frac{\partial}{\partial\tau^{2}}}\left( {{v_{3}^{s}\left( {\mu^{s},h^{s},\tau^{s}} \right)}{s\left( {t,\mu_{m\; i\; n}^{s},h^{s},\tau^{s}} \right)}} \right)} = {{- {\alpha^{s}\left( {\mu^{s},h^{s}} \right)}}{s\left( {t,\mu^{s},h^{s},\tau^{s}} \right)}}},} \\{{{{v^{s}\left( {E(t)} \right)}{s\left( {t,\mu_{m\; i\; n}^{s},h^{s},\tau^{s}} \right)}} = {r\left( {t,\mu_{{ma}\; x}^{r},h^{r},\tau^{r}} \right)}},} \\{{{s\left( {0,\mu^{s},h^{s},\tau^{s}} \right)} = {s_{0}\left( {\mu^{s},h^{s},\tau^{s}} \right)}},}\end{matrix} \right.$where s(t,μ^(s),h^(s),τ^(s)) is the population density of the marrowreticulocytes class at time t with maturity μ^(s), amount of hemoglobinh^(s) and number of TfRs τ^(s),13=μ_(min) ^(s)≦μ^(s)≦μ_(max) ^(s)=15, h_(min) ^(s)≦h^(s)≦h_(max) ^(s),τ_(min) ^(s)≦τ^(s)≦τ_(max) ^(s), t>0. The boundary conditionv^(s)(E(t))s(t,μ_(min) ^(s),h^(s),τ^(s))=r(t,μ_(max) ^(r),h^(r),τ^(r))describes the flux of cells leaving the erythroblast class and enteringthe reticulocytes population class carrying a certain amount ofhemoglobin and expressing a certain number of TfRs. Finally,s₀(μ^(s),h^(s),τ^(s)) is the population density at t=0.

A sigmoid function is used to describe the changes in the maturationvelocity v_(μ) ^(s),

${{v_{1}^{s}\left( {E(t)} \right)} = {\frac{a_{5} - b_{5}}{1 + {\mathbb{e}}^{{{- k_{5}}{E{(t)}}} + c_{5}}} + b_{5}}},$where E(t) is the EPO concentration at time t and a₅, b₅, c₅ and k₅ aresome positive constants with a₅>b₅. The maturation velocity increaseswith a rising concentration of EPO (Assumption 10). Note, a slowermaturation velocity causes the cells to reach the maximum cell ageμ_(max) ^(s) at a later point, whereas a faster maturation velocityshortens the transit time, i.e., the cells reach μ_(max) ^(s) earlier.The change of the amount of hemoglobin v_(h) ^(s) in a cell, thedevelopment of TfRs v_(τ) ^(s) on the cell surface and the rate ofapoptosis α^(s) are described with the same functions as for theerythroblast class, i.e.v ₂ ^(s)(τ,a ₁)=v ₂ ^(r)(τ,a ₁),v ₃ ^(s)(μ,h,τ)=v ₃^(r)(μ,h,τ),α^(s)(μ,h)=α^(r)(μ,h).

Further, β^(s)≡0 because of Assumption 22.

So far, Assumption 19 has not been accounted for. The assumption can beapplied in the following way: If E(t)>E* the reticulocyte class isskipped and erythroblasts with maximal maturity directly enter theerythrocytes class.

2.1.3. Erythrocytes.

Assumptions:

-   -   24. Erythrocytes and blood reticulocytes are subsumed in one        class.    -   25. Cells mature but do not proliferate.    -   26. There is a fixed random daily break-down of red blood cells        (not to be confused with loss of erythrocytes due to        senescence).    -   27. A drop in EPO concentration beneath a threshold precipitates        neocytolysis.    -   28. Erythrocytes with age between 14-21 days are likely to be        affected by neocytolysis.    -   29. Erythrocytes lose hemoglobin during senescence.    -   30. The expected life span of normochromic cells in healthy        persons is 120 days, but can decrease significantly in uremic        patients.    -   31. Cells in which hemoglobin content drops beneath a certain        threshold or which reach the maximum cell age are phagocytosed.

Up to this point, it has been assumed that red blood cells carryapproximately the same amount of hemoglobin. This assumption is notvalid in many pathologies, different types of anemia, and the majorityof dialysis patients, where a wide distribution of corpuscularhemoglobin can be observed. Additionally, it is desirable to account forhemoglobin which is lost during senescence. This loss, which amounts toabout 15-20%, is a non-negligible amount of iron that flows back to themacrophages. See S. C. Gifford, J. Derganc, Y. Shevkoplyas, S. S.Yoshida, and M. W. Bitensky, British Journal of Haematology,135:395-404, 2006; F. L. Willekens, F. H. Bosch, B.Roerdinkholder-Stoelwinder, Y. A. Groenen-Döpp, and J. M. Were, EuropeanJournal of Haematology, 58:246-250, 1997; J. M. Werre, F. Willekens, F.H. Bosch, L. D. de Haans, S. G. van der Vegt, A. G. van den Bos, and G.J. Bosman, Cellular and Molecular Biology 50 (2), 139-145, 2004; F. L.Willekens, B. Roerdinkholder-Stoelwinder, Y. A. Groenen-Döpp, H. J. Bos,G. J. Bosman, A. G. van den Bos, A. J. Verkleij, and J. M. Were, Blood,101:747-751, 2003. Further, the binding of oxygen is pivotal fordelivery of oxygen by hemoglobin, and thus a reduction of 15-20% in thehemoglobin mass will have a noticeable impact on the oxygen carryingcapacity. Hence, there are several good reasons that suggest includinghemoglobin content h^(m) as a structural variable. Further, theattribute h^(m) allows one to relatively easily compute the amount ofiron which is set free when an erythrocyte is phagocytosed. Moreover, itenables tracking of how much iron is set free daily because of sheddingof hemoglobin containing vesicles.

Altogether, the following equation is obtained for erythrocytes:

$\left\{ \begin{matrix}{{\left. {\frac{\partial\;}{{\partial t}\;}{m\left( {t,\mu^{m},h^{m}} \right)}} \right) + {\frac{\partial\;}{\partial\mu^{m}}{m\left( {t,\mu^{m},h^{m}} \right)}} + {{v_{2}^{m}\left( \mu^{m} \right)}\frac{\partial\;}{\partial h^{m}}{m\left( {t,\mu^{m},h^{m}} \right)}}},} \\{{= {- {\alpha^{m}\left( {{E(t)},\mu^{m},h^{m}} \right)}}},{m\left( {t,\mu^{m}} \right)},} \\{{{m\left( {t,\mu^{m},h_{\min}^{m}} \right)} = 0},} \\{{m\left( {t,0,h^{m}} \right)} = \left\{ \begin{matrix}{{v^{s}\left( {E(t)} \right)}{\int_{\tau_{\min}^{s}}^{\tau_{\max}^{s}}{{s\left( {t,\mu_{\max}^{s},h^{s},\tau^{s}} \right)}\ {\mathbb{d}\tau^{s}}}}} & {{{{for}\mspace{14mu}{E(t)}} \leq E^{*}},} \\{\int_{\tau_{\min}^{r}}^{\tau_{\max}^{r}}{{r\left( {t,\mu_{\max}^{r},h^{r},\tau^{r}} \right)}\ {\mathbb{d}\tau^{r}}}} & {{{{for}\mspace{14mu}{E(t)}} > E^{*}},}\end{matrix}\quad \right.} \\{{{m\left( {0,\mu^{m},h^{m}} \right)} = {m_{0}\left( {\mu^{m},h^{m}} \right)}},}\end{matrix}\quad \right.$where m(t,μ^(m),h^(m)) is the population density for the erythrocyteclass at time t with maturity μ^(m) and cell hemoglobin h^(m), 0=μ_(min)^(m)≦μ^(m)≦μ_(max) ^(m)=120, h_(min) ^(m)≦h^(m)≦h_(max) ^(m) and t≧0 andm₀(μ^(m),h^(m)) is the population density at t=0. Moreover, m(t,0,h^(m))is defined by the number of reticulocytes or the number of erythroblastsentering the blood stream carrying a certain amount of hemoglobin h^(m),respectively. To describe the corpuscular hemoglobin decay, a functionv_(h) ^(m) is used, of the formv ₂ ^(m)(μ^(m))=b ₆(μ^(m)−μ_(*) ^(m))² −c ₆,with b₆, c₆, μ_(*) ^(m)>0 constant. Choosing appropriate constants b₆,c₆ and μ_(*) ^(m) allows for a wide range of possible forms of thefunction h^(m)(t), with dh^(m)(t)/dt=v₂ ^(m)(μ^(m)(t)).

The mortality rate α^(m)(E(t),μ^(m),h^(m)) is the sum of three terms:

-   -   1. The rate of cells dying due to senescence, which is assumed        to happen when the cell reaches the maximal cell age μ_(max)        ^(m) or the minimal hemoglobin content h_(min) ^(m) (Assumption        30).    -   2. Cells dying because of random break-down.    -   3. Cells dying because of neocytolysis.        Altogether,

${\alpha^{m}\left( {{E(t)},\mu^{m},h^{m}} \right)} = \left\{ \begin{matrix}{{\gamma^{m}\left( {\mu^{m},h^{m}} \right)} + {\min\left( {{c_{E}/{E(t)}^{k_{E}}},b_{E}} \right)}} & {{{{for}\mspace{14mu}{E(t)}} < \tau_{E}},{\mu_{\min}^{m,n} \leq \mu^{m} \leq \mu_{\max}^{m,n}},} \\{\gamma^{m}\left( {\mu^{m},h^{m}} \right)} & {{otherwise},}\end{matrix} \right.$where b_(E), c_(E), k_(E) are positive constants, τ_(E) is the thresholdbeneath which neocytolysis is triggered (Assumption 27). Further,[μ_(min) ^(m,n),μ_(max) ^(m,n)]=[14,21] is the age-interval during whichcells are affected by neocytolysis (Assumption 28) and t≧0. Since it isassumed that a cell is phagocytosed when it reaches the maximal cell ageμ_(max) ^(m) or the minimal cell hemoglobin h_(min) ^(m), the mortalityrate γ^(m)(μ^(m),h^(m)) needs to be chosen such thatγ^(m)(μ^(m),h^(m))=α_(rand) ^(m) for μ_(min) ^(m)≦μ^(m)≦μ_(max) ^(m)−δ,h_(min) ^(m)+δ≦h^(m)≦h_(max) ^(m) with δ>0 sufficiently small and ∫_(μ)_(max) _(−δ) ^(μ) ^(max) ∫_(h) _(min) ^(h) ^(min) ^(+δ)α^(m)(μ,h)dμdh=∞.Here, α_(rand) ^(m) denotes the random breakdown of cells (Assumption26). A possible choice for γ^(m)(μ^(m),h^(m)) is

${\gamma^{m}\left( {\mu^{m},h^{m}} \right)} = \left\{ {{{\begin{matrix}\alpha_{r}^{m} & {{{{for}\mspace{14mu}\mu^{m}} \in \left\lbrack {\mu_{\min}^{m},{\mu_{\max}^{m} - \delta}} \right\rbrack},{h^{m} \in \left\lbrack {{h_{\min}^{m} + \delta},h_{\max}^{m}} \right\rbrack},} \\{{\gamma_{\mu}^{m}\left( \mu^{m} \right)}{\gamma_{h}^{m}\left( h^{m} \right)}} & {{{{for}\mspace{14mu}\mu^{m}} \in \left( {{\mu_{\max}^{m} - \delta},\mu_{\max}^{m}} \right)},{h^{m} \in \left( {h_{\min}^{m},{h_{\min}^{m} + \delta}} \right)},} \\\infty & {{{{for}\mspace{14mu}\mu^{m}} \geq \mu_{\max}^{m}},{h^{m} \leq h_{\min}^{m}},}\end{matrix}\mspace{79mu}{where}\mspace{79mu}{\gamma_{\mu}^{m}\left( \mu^{m} \right)}} = \frac{3\alpha_{rand}^{m}\delta^{2}}{\left( \mu^{m} \right)^{2} - {2\left( {\mu_{\max}^{m} + \delta} \right)\mu^{m}} + {\left( {\mu_{\max}^{m} + {2\delta}} \right)\mu_{\max}^{m}}}},\mspace{79mu}{{{and}\mspace{79mu}{\gamma_{h}^{m}\left( h^{m} \right)}} = \frac{3\delta^{2}}{\left( h^{m} \right)^{2} - {2\left( {h_{\min}^{m} - \delta} \right)h^{m}} + {\left( {h_{\min}^{m} - {2\delta}} \right)h_{\max}^{m}}}},} \right.$In case of bleeding an additional term α_(bleed) ^(m)(t) needs to beadded. What this term looks like depends on whether the bleeding lastsminutes or hours or is prolonged for several days or weeks.2.1.4. Erythropoietin.Assumptions:

-   -   32. Release of EPO is controlled by a negative feedback        mechanism according to the oxygen content.    -   33. Oxygen carrying capacity is directly proportional to the        amount of hemoglobin of circulating red blood cells.    -   34. The degradation rate of EPO is constant.    -   35. There is a slight delay in reaction of the EPO production        rate to the number of RBCs but this is negligible compared to        the duration of development of erythrocytes.

The equation which describes the release of EPO by the kidneys is

${{E_{in}^{end}(t)} = {\frac{a_{7} - b_{7}}{1 + {\mathbb{e}}^{{k_{7}\overset{\sim}{H}{b{(t)}}} - c_{7}}} + b_{7}}},$where {tilde over (H)}b(t)=Hb(t)/TBV is the hemoglobin concentration.Here, TBV is the total blood volume and Hb(t)=∫_(h) _(min) _(m) ^(h)^(max) ^(m) m(t,μ^(m),h^(m))dh^(m) is the total amount of hemoglobin inblood. Furthermore, a₇, b₇, c₇ and k₇ are positive constants with a₇>b₇.

The dynamics of the endogenous EPO concentration E^(end)(t) aredescribed by the following ordinary differential equation:

${{{\overset{.}{E}}^{end}(t)} = {{\frac{1}{TBV}{E_{in}^{{end}\;}(t)}} - {c_{\deg}^{end}{E^{end}(t)}}}},$where E^(end)(t) is the endogenous EPO concentration, E_(in) ^(end) isthe amount of EPO released by the kidneys and c_(deg) ^(end) describesthe degradation rate of endogenous EPO. The change in the concentrationE^(ex)(t) of an ESA reads

${{{\overset{.}{E}}^{ex}(t)} = {{\frac{1}{TBV}{E_{in}^{{ex}\;}(t)}} - {c_{\deg}^{ex}{E^{ex}(t)}}}},$where E_(in) ^(ex)(t) is the rate at which the artificial hormone isadministered c_(deg) ^(ex) is the rate with which the exogenous hormoneis degraded. In intravenous administration, the total amount of theagent is injected into a vein, within a very short time interval. Inthis case E_(in) ^(ex)(t) can be approximated by E₀ ^(ex)(t)δ_(t) ₀ (t),where E₀ ^(ex) is the amount of artificial hormone administered andδ_(t) ₀ (t) is the Dirac delta impulse located at t₀, the time when theadministration takes place. In addition, the degradation rate forexogenous EPO c_(deg) ^(ex) differs from the one for endogenous EPO andvaries according to the kind of ESA administered. The overallconcentration of EPO in plasma consists of the naturally producederythropoietin in the body and the administered ESAE(t)=E ^(ex)(t)+E ^(end)(t).3. The Iron Submodel.Assumptions:

-   -   37. Five iron pools are considered within the body: a plasma        pool, iron in precursor cells, iron in erythrocytes, iron in the        macrophages and a storage pool.    -   38. Other iron in the body is ignored.    -   39. Hepcidin regulates iron homeostasis.    -   40. Stressed erythropoiesis and depleted stores decrease the        amount of hepcidin released into plasma.    -   41. Large iron stores and inflammation increase hepcidin levels.    -   42. The amount of iron in plasma is proportional to the amount        of transferrin molecules carrying iron. (Free iron is        neglected.)    -   43. No distinction is made between monoferric and differic        transferrin.    -   44. Iron is lost via loss of cells.    -   45. Iron lost via urine and sweat is neglected.

Iron is a part of hemoglobin which provides the means of O₂ transport tothe tissue. The rate at which iron can be released to plasma fromexisting stores and macrophages may easily limit the rate of irondelivery for hemoglobin synthesis. Therefore, a strong relationshipbetween impaired erythropoiesis and total body iron burden exists. Anexcellent description of iron homeostasis and its effects can be foundin Crichton.

Iron is a substance which is very tightly controlled from the body asiron overload is toxic. Iron homeostasis can be only achieved by thecontrol of absorption because the human body, unlike other mammals, isnot able to influence the excretion of iron.

Under normal circumstances the average daily loss of iron is very small(men: 1 mg, women: 2 mg, because of menstruation or pregnancy).

For the iron submodel, five dynamic iron pools are considered within thebody: a plasma pool, iron in precursor cells, iron in erythrocytes, ironin the macrophages and a storage pool. Hence, ferritin and hemosiderinstores are combined. Other iron in the body is ignored. For anillustration of the organizational structure of the iron submodel seeFIG. 8. The amount of iron lost via urine and sweat is neglected. Thus,in this model, excretion of iron takes place only when cells are lost.On one hand, this can be due to loss of RBCs, i.e., due to externalbleeding, and on the other hand there is a daily loss of epithelialcells, which amount for about 1 mg/day. Further, it is assumed that onlythe precursor cells in the bone marrow uptake iron and ineffectiveerythropoiesis depends on the amount of iron which can be provided forerythropoiesis. Moreover, the model accounts for the loss of iron fromerythrocytes during senescence.

The general form of the compartmental iron model is as follows:

${{\frac{\mathbb{d}}{\mathbb{d}t}a_{1}} = {{{k_{41}(H)}a_{4}} + {k_{{iv},1}a_{iv}} + {{k_{{gastro},1}(H)}a_{gastro}} + {{k_{51}(H)}a_{5}} - {k_{15}a_{1}} - {k_{12}a_{1}}}},\mspace{79mu}{{\frac{\mathbb{d}}{\mathbb{d}t}a_{2}} = {{k_{12}a_{1}} - f_{24} - f_{23}}},\mspace{79mu}{{\frac{\mathbb{d}}{\mathbb{d}t}a_{3}} = {f_{23} - f_{34} - f_{3,{out}}}},\mspace{79mu}{{\frac{\mathbb{d}}{\mathbb{d}t}a_{4}} = {f_{24} + {f_{34}a_{3}} - {{k_{41}(H)}a_{4}} - {k_{45}a_{4}}}},\mspace{79mu}{{\frac{\mathbb{d}}{\mathbb{d}t}a_{5}} = {{k_{15}a_{1}} + {k_{45}a_{4}} - {{k_{51}(H)}{a_{5}.}}}}$Note, in the above equations, the dependence of particular components ont was not explicitly stated to simplify notation. Furthermore, a_(i),i=1, . . . , 5, denotes the amount of iron in the compartment i andk_(ij), i, j=1, . . . , 5 are the transfer rates from compartment i tocompartment j, i.e., k_(ij)a_(i) is the rate at which iron is movingfrom compartment i to compartment j. Further, ƒ_(ij) denote the rates atwhich iron is transferred from compartment i into compartment j, when itis not of the simple form k_(ij)a_(i); ƒ_(3,out) is the amount of ironwhich is lost by bleeding. Additionally, a_(iv) denotes the amount ofiron intravenously administered, whereas a_(gastro) is the amount ofiron in the duodenum. Finally, k_(iv,1) and k_(gastro,1)(H) are thecorresponding flow rates. The rate k_(gasto,1) can be described by adecreasing sigmoidal function

${{k_{{gastro},1}\left( {H(t)} \right)} = {\frac{a_{8} - b_{8}}{1 + {\mathbb{e}}^{{{- k_{8}}{H{(t)}}} + c_{8}}} + b_{8}}},$where a₈, b₈, c₈ and k₈ are positive constants with a₈>b₈. The functionH(t) is the solution of the differential equation (5).

Since the time needed to administer iron is very short compared to othertime constants in the system, the following can be takenk _(iv,1) a _(iv) =a _(iv,total)δ_(t) ₀ .Here a_(iv,total) is the total amount of iron administered and δ_(t) ₀is the Dirac delta function located at t₀.3.1. Hepcidin Feedback.

The regulator of iron homeostasis is hepcidin, a peptide hormoneproduced by the liver. Hepcidin negatively regulates the availability ofiron in plasma by binding to ferroportin on the membranes ofiron-exporting cells. Hepcidin binding induces the endocytosis andproteolysis of ferroportin, and thereby decreases the delivery of ironto the plasma. See Nemeth E., Tuttle M. S., Powelson J., Vaughn M. B.,Donovan A., Ward D. M., Ganz T., and Kaplan J., Science 306 pp.2090-2093 (2004). Ferroportin is the iron exporter required for ironegress from iron exporting cells, as for instance, enterocytes andmacrophages. Hepcidin acts as a medium to withhold iron from certaininvasive bacteria, as they are not able to proliferate under conditionsof insufficient iron supply. The dynamics of the hepcidin concentrationin plasma is calculated according to

$\begin{matrix}{{{\frac{\mathbb{d}\;}{\mathbb{d}t}{H(t)}} = {{\frac{1}{TBV}{H_{prod}(t)}} - {c_{\deg}^{H}{H(t)}}}},} & (5)\end{matrix}$where H_(prod)(t) is the rate at which hepcidin, produced by thehepatocytes, is released into the plasma at time t. TBV is the totalblood volume, and c_(deg) ^(H) is the rate at which hepcidin isdegraded. The production rate H_(prod)(t) of hepcidin is obtained as theresult of a feedback loop involving the total precursor cell population,the amount of iron stored and the state of inflammation of the patientH _(prod)(t)=ƒ(U(t),a ₅(t),i(t)).Here

${{U(t)} = {{\int\limits_{\mu_{\min}^{r}}^{\mu_{\max}^{r}}{\int\limits_{h_{\min}^{r}}^{h_{\max}^{r}}{\int\limits_{\tau_{\min}^{r}}^{\tau_{\max}^{r}}{{r\left( {t,\mu^{r},h^{r}} \right)}{\mathbb{d}\mu^{r}}{\mathbb{d}h^{r}}{\mathbb{d}\tau^{r}}}}}} + {\int\limits_{\mu_{\min}^{s}}^{\mu_{\max}^{s}}{\int\limits_{h_{\min}^{s}}^{h_{\max}^{s}}{\int\limits_{\tau_{\min}^{s}}^{\tau_{\max}^{s}}{{s\left( {t,\mu^{s},h^{s}} \right)}{\mathbb{d}\mu^{s}}{\mathbb{d}h^{s}}{\mathbb{d}\tau^{s}}}}}}}},$where U(t) describes the total precursor cell population(erythroblasts+reticulocytes). Further, a₅(t) is the amount of ironstored and i(t) is the inflammation status at time t. The feedback ismodeled using a monotonically increasing sigmoidal function

${{H_{prod}(t)} = {\frac{a_{H} - b_{H}}{1 + {\mathbb{e}}^{{{- k_{H}}{C{(t)}}} + c_{H}}} + b_{H}}},$where a_(H), b_(H), c_(H) and k_(H) are positive constants witha_(H)>b_(H). The functionC(t)=−c _(U) U(t)+c _(a) ₅ a ₅(t)+c _(i) i(t),is an auxiliary equation that increases if storage amount orinflammation status rises, and decreases when more RBCs are produced.Here, c_(U), c_(a) ₅ and c_(i) are positive constants.3.2. Detailed Description of the Iron Compartments.3.2.1. Plasma (Compartment 1).

The rate at which iron enters erythroid cells at time t is given by

${d_{0}{a_{1}(t)}\left( {{\int\limits_{\mu_{\min}^{r}}^{\mu_{\max}^{r}}{\int\limits_{h_{\min}^{r}}^{h_{\max}^{r}}{\int\limits_{\tau_{\min}^{r}}^{\tau_{\max}^{r}}{\tau^{r}{r\left( {t,\mu^{r},h^{r}} \right)}{\mathbb{d}\mu^{r}}{\mathbb{d}h^{r}}{\mathbb{d}\tau^{r}}}}}} + {\int\limits_{\mu_{\min}^{s}}^{\mu_{\max}^{s}}{\int\limits_{h_{\min}^{s}}^{h_{\max}^{s}}{\int\limits_{\tau_{\min}^{s}}^{\tau_{\max}^{s}}{\tau^{s}{s\left( {t,\mu^{s},h^{s}} \right)}{\mathbb{d}\mu^{s}}{\mathbb{d}h^{s}}{\mathbb{d}\tau^{s}}}}}}} \right)},$where r(t,μ^(r),τ^(r),h^(r)) and s(t,μ^(s),τ^(s),h^(s)) denote thepopulation densities of hemoglobin synthesizing cells, erythroblasts andreticulocytes, respectively. At the same time this is the rate at whichiron leaves the plasma compartment and enters the precursor cellcompartment.

Only a very small amount of iron can be found in this pool (about 3 mg;≦1%) bound to transferrin, but it is a very important compartmentbecause of its rapid turnover rate. See Finch 1982. Iron normally turnsover at least 10 times each day. See Williams Hematology. Inflow intothe plasma compartment comes from macrophages, the storage pool, ironabsorbed through the duodenum and via intravenously administered iron.Outflow leaves the plasma compartment to the precursor cells and thestorage compartment. The flow rates k₄₁, k_(gastro,1) and k₅₁ aredependent on the hepcidin concentration H(t) at time t. They decreasewhen the hepcidin level rises and vice versa. Inflammation and full ironstores increase the release of hepcidin from the hepatocytes in theliver, whereas a stressed erythropoiesis and depleted stores suppressthe production of hepcidin. The outflow to the precursor cellcompartment k₁₂a₁ is altered by the amount of iron needed from the cellsto synthesize hemoglobin. As erythroid cells are only able to take upiron through transferrin receptors, their need for iron can beassociated with the number of receptors expressed on precursor cells.Thus, the rate k₁₂ is given by

${k_{12} = {d_{0}\left( {{\int\limits_{\mu_{\min}^{r}}^{\mu_{\max}^{r}}{\int\limits_{h_{\min}^{r}}^{h_{\max}^{r}}{\int\limits_{\tau_{\min}^{r}}^{\tau_{\max}^{r}}{\tau^{r}{r\left( {t,\mu^{r},h^{r}} \right)}{\mathbb{d}\mu^{r}}{\mathbb{d}h^{r}}{\mathbb{d}\tau^{r}}}}}} + {\int\limits_{\mu_{\min}^{s}}^{\mu_{\max}^{s}}{\int\limits_{h_{\min}^{s}}^{h_{\max}^{s}}{\int\limits_{\tau_{\min}^{s}}^{\tau_{\max}^{s}}{\tau^{s}{s\left( {t,\mu^{s},h^{s}} \right)}{\mathbb{d}\mu^{s}}{\mathbb{d}h^{s}}{\mathbb{d}\tau^{s}}}}}}} \right)}},$For an explanation of the various terms see the statements concerningerythroblasts in Subsection 2.1.2. Hence, the change of amount of ironin the plasma compartment is

${\frac{\mathbb{d}\;}{\mathbb{d}t}{a_{1}(t)}} = {{{k_{41}\left( {H(t)} \right)}{a_{4}(t)}} + {{a_{{iv},{total}}(t)}{\delta_{t_{0}}(t)}} + {{k_{{gastro},1}\left( {H(t)} \right)}{a_{gastro}(t)}} + {{k_{51}\left( {H(t)} \right)}{a_{5}(t)}} - {k_{15}{a_{1}(t)}} - {{d_{0}\left( {{\int\limits_{\mu_{\min}^{r}}^{\mu_{\max}^{r}}{\int\limits_{h_{\min}^{r}}^{h_{\max}^{r}}{\int\limits_{\tau_{\min}^{r}}^{\tau_{\max}^{r}}{\tau^{r}{r\left( {t,\mu^{r},h^{r}} \right)}{\mathbb{d}\mu^{r}}{\mathbb{d}h^{r}}{\mathbb{d}\tau^{r}}}}}} + {\int\limits_{\mu_{\min}^{s}}^{\mu_{\max}^{s}}{\int\limits_{h_{\min}^{s}}^{h_{\max}^{s}}{\int\limits_{\tau_{\min}^{s}}^{\tau_{\max}^{s}}{\tau^{s}{s\left( {t,\mu^{s},h^{s}} \right)}{\mathbb{d}\mu^{s}}{\mathbb{d}h^{s}}{\mathbb{d}\tau^{s}}}}}}} \right)}{{a_{1}(t)}.}}}$The transfer rate k₄₁(H(t)) and k₅₁(H(t)) can be described bymonotonically decreasing sigmoidal functions

$\begin{matrix}{{{k_{41}\left( {H(t)} \right)} = {\frac{a_{9} - b_{9}}{1 + {\mathbb{e}}^{{{- k_{9}}{H{(t)}}} + c_{9}}} + b_{9}}},{and}} & (6) \\{{{k_{51}\left( {H(t)} \right)} = {\frac{a_{10} - b_{10}}{1 + {\mathbb{e}}^{{{- k_{10}}{H{(t)}}} + c_{10}}} + b_{10}}},} & (7)\end{matrix}$respectively. The parameters a₉, b₉, c₉, k₉, a₁₀, b₁₀, c₁₀ and k₁₀ arepositive constants with a₉>b₉ and a₁₀>b₁₀. The hepcidin level H(t) attime t is determined by equation (5).3.2.2. Precursor Cells (Compartment 2).

Inflow into the precursor cell compartment comes from plasma. Oneoutflow of this compartment leaves to the erythrocytes compartment anddescribes the amount of hemoglobin carried from reticulocytes, whichenter the blood stream. Thus, the flow ƒ₂₃(•) from compartment 2 tocompartment 3, is defined by the cells that leave the bone marrow andthe amount of hemoglobin these cells are carrying. These quantities canbe calculated from the population model andk ₂₃ a ₂=ƒ₂₃(v ^(s)(E(t))s(t,μ _(max) ^(s),τ^(s)),h ^(s)(μ_(max) ^(s))),where v^(s)(E(t))s(t,μ_(max) ^(s),τ^(s)) determines how many cells areleaving the bone marrow at a certain time t and the distributionfunction h^(s)(μ_(max) ^(s)) at the maximum maturity level μ_(max) ^(s)indicates how much hemoglobin they carry. Here v^(s)(E(t)) is thematuration velocity of reticulocytes depending on EPO, s(t,μ_(max)^(s),τ^(s)) is the density of reticulocyte population with maximalmaturity μ_(max) ^(s), and τ^(s)ε[τ_(min) ^(s),τ_(max) ^(s)] is thenumber of transferrin receptors at time t. The hemoglobin distributionfunction h^(s)(μ_(max) ^(s)) arises as a result of the population model,in particular the precursor cell population part.

The second outflow of this compartment leaves to the macrophages and isdefined by the precursor cells that die. Hence, the flow fromcompartment 2 to compartment 4 is defined by the amount of iron thatdying cells carry, i.e.,

$\begin{matrix}{{{k_{24}a_{4}} = {{\int_{\mu_{\min}^{r}}^{\mu_{\max}^{r}}{\int_{h_{\min}^{r}}^{h_{\max}^{r}}{\int_{\tau_{\min}^{r}}^{\tau_{\max}^{r}}{h^{r}{\alpha^{r}\left( {\mu^{r},h^{r}} \right)}{r\left( {t,\mu^{r},h^{r}} \right)}\ {\mathbb{d}\mu^{r}}{\mathbb{d}h^{r}}{\mathbb{d}\tau^{r}}}}}} + {\int_{\mu_{\min}^{s}}^{\mu_{\max}^{s}}{\int_{h_{\min}^{s}}^{h_{\max}^{s}}{\int_{\tau_{\min}^{s}}^{\tau_{\max}^{s}}{h^{s}{\alpha^{s}\left( {\mu^{s},h^{s}} \right)}{s\left( {t,\mu^{s},h^{s}} \right)}\ {\mathbb{d}\mu^{s}}{\mathbb{d}h^{s}}{\mathbb{d}\tau^{s}}}}}}}},} & (8)\end{matrix}$where r(t,μ^(r),τ^(r)) and s(t,μ^(s),h^(s)) are the densities of the twoprecursor cell populations. Further, μ^(i)ε[μ_(min) ^(i),μ_(max) ^(i)]denotes the maturity, h^(i)ε[h_(min) ^(i),h_(max) ^(i)] denotes theamount of hemoglobin carried by a cell and τ^(i)ε[τ_(min) ^(i),τ_(max)^(i)] denotes the number of transferrin receptors for i=r, s. Thefunction α^(i)(μ^(i),h^(i)), i=r, s describes the mortality rate ofprecursor cells, where α^(r)=α^(s) is determined by equations (3) and(4).

Altogether the change of amount of iron in the iron-precursorscompartment is given by

${\frac{\mathbb{d}}{\mathbb{d}t}{a_{2}(t)}} = {{{d_{0}\left( {{\int_{\mu_{\min}^{r}}^{\mu_{\max}^{r}}{\int_{h_{\min}^{r}}^{h_{\max}^{r}}{\int_{\tau_{\min}^{r}}^{\tau_{\max}^{r}}{\tau^{r}{r\left( {t,\mu^{r},h^{r}} \right)}\ {\mathbb{d}\mu^{r}}{\mathbb{d}h^{r}}{\mathbb{d}\tau^{r}}}}}} + {\int_{\mu_{\min}^{s}}^{\mu_{\max}^{s}}{\int_{h_{\min}^{s}}^{h_{\max}^{s}}{\int_{\tau_{\min}^{s}}^{\tau_{\max}^{s}}{\tau^{s}{s\left( {t,\mu^{s},h^{s}} \right)}\ {\mathbb{d}\mu^{s}}{\mathbb{d}h^{s}}{\mathbb{d}\tau^{s}}}}}}} \right)}{a_{1}(t)}} - {f_{23}\left( {{{v^{s}\left( {E(t)} \right)}{s\left( {t,\mu_{\max}^{s},\tau^{s}} \right)}},{h^{s}\left( \mu_{\max}^{s} \right)}} \right)} - {\int_{\mu_{\min}^{r}}^{\mu_{\max}^{r}}{\int_{h_{\min}^{r}}^{h_{\max}^{r}}{\int_{\tau_{\min}^{r}}^{\tau_{\max}^{r}}{h^{r}{\alpha^{r}\left( {\mu^{r},h^{r}} \right)}{r\left( {t,\mu^{r},h^{r}} \right)}\ {\mathbb{d}\mu^{r}}{\mathbb{d}h^{r}}{\mathbb{d}\tau^{r}}}}}} + {\int_{\mu_{\min}^{s}}^{\mu_{\max}^{s}}{\int_{h_{\min}^{s}}^{h_{\max}^{s}}{\int_{\tau_{\min}^{s}}^{\tau_{\max}^{s}}{h^{s}{\alpha^{s}\left( {\mu^{s},h^{s}} \right)}{s\left( {t,\mu^{s},h^{s}} \right)}\ {\mathbb{d}\mu^{s}}{\mathbb{d}h^{s}}{{\mathbb{d}\tau^{s}}.}}}}}}$Note that it is not necessary to determine the function ƒ₂₃, becausethere is no need to explicitly observe the rate of change of theprecursor cell compartment. The amount of iron a₂(t) in this compartmentat time t can be directly computed from the population classes. It isdefined by the total number of precursor cells and the amount ofhemoglobin they carry:

a₂(t) = ∫_(μ_(min)^(r))^(μ_(max)^(r))∫_(h_(min)^(r))^(h_(max)^(r))∫_(τ_(min)^(r))^(τ_(max)^(r))h^(r)r(t, μ^(r), h^(r)) 𝕕μ^(r)𝕕h^(r)𝕕τ^(r) + ∫_(μ_(min)^(s))^(μ_(max)^(s))∫_(h_(min)^(s))^(h_(max)^(s))∫_(τ_(min)^(s))^(τ_(max)^(s))h^(s)s(t, μ^(s), h^(s)) 𝕕μ^(s)𝕕h^(s)𝕕τ^(s).

However, it is important to know the flux of hemoglobin to themacrophage compartment which arises from dying cells, because this flux(see equation (8)) is a source term in the macrophage compartment.

3.2.3. Erythrocytes (Compartment 3).

Inflow into the erythrocytes compartment comes from the precursor cellcompartment, and iron leaves the compartment to the macrophagescompartment (phagocytosis of senescent cells, random break-down of cellsand neocytolysis, shedding of hemoglobin containing vesicles)

$\begin{matrix}{{{k_{34}a_{3}} = {{\int_{\mu_{\min}^{m}}^{\mu_{\max}^{m}}{\int_{h_{\min}^{m}}^{h_{\max}^{m}}{h^{m}{\alpha^{m}\left( {{E(t)},\mu^{m},h^{m}} \right)}{m\left( {t,\mu^{m},h^{m}} \right)}\ {\mathbb{d}\mu^{m}}\ {\mathbb{d}h^{m}}}}} + {\int_{\mu_{\min}^{m}}^{\mu_{\max}^{m}}{\int_{h_{\min}^{m}}^{h_{\max}^{m}}{{v_{2}^{m}\left( \mu^{m} \right)}{m\left( {t,\mu^{m},h^{m}} \right)}\ {\mathbb{d}\mu^{m}}\ {\mathbb{d}h^{m}}}}}}},} & (9)\end{matrix}$or is excreted via bleeding

k_(3, out)a₃ = ∫_(μ_(min)^(m))^(μ_(max)^(m))∫_(h_(min)^(m))^(h_(max)^(m))h^(m)α_(bleed)^(m)(t)m(t, μ^(m), h^(m)) 𝕕μ^(m) 𝕕h^(m).Here α^(m) denotes the mortality rate of erythrocytes due tophagocytosis of senescent cells, random break-down of cells andneocytolysis. Furthermore, h^(m) is the amount of hemoglobin carried bya cell, m(t,μ^(m),h^(m)) is the population density of erythrocytes withcell age μ^(m) and cellular hemoglobin h^(m) at time t and α_(bleed)^(m)(t) is the rate with which RBCs are lost due to bleeding. Thus, thechange of amount of iron in the erythrocyte compartment is

${\frac{\mathbb{d}}{\mathbb{d}t}{a_{3}(t)}} = {{f_{23}\left( {{{v^{s}\left( {E(t)} \right)}{s\left( {t,\mu_{\max}^{s},\tau^{s}} \right)}},{h^{s}\left( \mu_{\max}^{s} \right)}} \right)} - {\int_{\mu_{\min}^{m}}^{\mu_{\max}^{m}}{\int_{h_{\min}^{m}}^{h_{\max}^{m}}{h^{m}{\alpha_{bleed}^{m}(t)}{m\left( {t,\mu^{m},h^{m}} \right)}\ {\mathbb{d}\mu^{m}}\ {\mathbb{d}h^{m}}}}} - {\int_{\mu_{\min}^{m}}^{\mu_{\max}^{m}}{\int_{h_{\min}^{m}}^{h_{\max}^{m}}{h^{m}{\alpha^{m}\left( {{E(t)},\mu^{m},h^{m}} \right)}{m\left( {t,\mu^{m},h^{m}} \right)}\ {\mathbb{d}\mu^{m}}\ {\mathbb{d}h^{m}}}}} - {\int_{\mu_{\min}^{m}}^{\mu_{\max}^{m}}{\int_{h_{\min}^{m}}^{h_{\max}^{m}}{{v_{2}^{m}\left( \mu^{m} \right)}{m\left( {t,\mu^{m},h^{m}} \right)}\ {\mathbb{d}\mu^{m}}\ {{\mathbb{d}h^{m}}.}}}}}$Similar to compartment 2, it is not necessary to observe the rate ofchange of the erythrocyte cell compartment to compute the amount of irona₃(t) at time t. The quantity a₃(t) can be directly computed from thepopulation class and is defined by the total number of circulating redblood cells and the amount of hemoglobin they carry:

a₃(t) = ∫_(μ_(min)^(m))^(μ_(max)^(m))∫_(h_(min)^(m))^(h_(max)^(m))h^(m)m(t, μ^(m), h^(m)) 𝕕h^(m) 𝕕μ^(m).However, the term k₃₄a₃(t) needs to be known explicitly (see equation(9)). This flux, which is the rate of hemoglobin released byerythrocytes, appears as a source term in the macrophage compartment.3.2.4. Macrophages (Compartment 4).

Inflow into the macrophages compartment comes from the erythrocytescompartment (erythrocytes phagocytosed by macrophages, iron set free byrandom break-down of RBCs and shedding of hemoglobin containingvesicles) and the precursor cells compartment. Outflow of thiscompartment leaves to the plasma and to the storage compartment.Therefore, the change of iron in the compartment is

${{\frac{\mathbb{d}}{\mathbb{d}t}{a_{4}(t)}} = {{\int_{\mu_{\min}^{m}}^{\mu_{\max}^{m}}{\int_{h_{\min}^{m}}^{h_{\max}^{m}}{h^{m}{\alpha^{m}\left( {{E(t)},\mu^{m},h^{m}} \right)}{m\left( {t,\mu^{m},h^{m}} \right)}\ {\mathbb{d}\mu^{m}}\ {\mathbb{d}h^{m}}}}} + {\int_{\mu_{\min}^{m}}^{\mu_{\max}^{m}}{\int_{h_{\min}^{m}}^{h_{\max}^{m}}{{v_{2}^{m}\left( \mu^{m} \right)}h^{m}{m\left( {t,\mu^{m},h^{m}} \right)}\ {\mathbb{d}\mu^{m}}\ {\mathbb{d}h^{m}}}}} + {\int_{\mu_{\min}^{r}}^{\mu_{\max}^{r}}{\int_{h_{\min}^{r}}^{h_{\max}^{r}}{\int_{\tau_{\min}^{r}}^{\tau_{\max}^{r}}{h^{r}{r\left( {t,\mu^{r},h^{r}} \right)}\ {\mathbb{d}\mu^{r}}{\mathbb{d}h^{r}}{\mathbb{d}\tau^{r}}}}}} + {\int_{\mu_{\min}^{s}}^{\mu_{\max}^{s}}{\int_{h_{\min}^{s}}^{h_{\max}^{s}}{\int_{\tau_{\min}^{s}}^{\tau_{\max}^{s}}{h^{s}{s\left( {t,\mu^{s},h^{s}} \right)}\ {\mathbb{d}\mu^{s}}{\mathbb{d}h^{s}}{\mathbb{d}\tau^{s}}}}}} - {{k_{41}\left( {H(t)} \right)}{a_{4}(t)}} - {k_{45}{a_{4}(t)}}}},$where k₄₁(H(t)) is given by equation (6).3.2.5. Storage (Compartment 5.)

Inflow into the storage compartment comes from the plasma and themacrophages compartment and iron leaves this compartment to the plasmacompartment. Hence, the change of iron in the storage compartment is

${{\frac{\mathbb{d}}{\mathbb{d}t}{a_{5}(t)}} = {{k_{15}{a_{1}(t)}} + {k_{45}{a_{4}(t)}} - {{k_{51}\left( {H(t)} \right)}{a_{5}(t)}}}},$where k₅₁(H(t)) is given in equation (7).

The iron submodel can be applied independently of erythropoiesis to avariety of iron disorders, such as true iron deficiency states (e.g.,bleeding, malabsorption), functional iron deficiencies (e.g., anemias ofchronic disease, usually with chronic inflammatory disorders ormalignancies), or iron overload disorders (e.g., all forms ofhemochromatosis, iron intoxication, polytransfusions of blood, chronichepatopathies, and hyper plastic erythroid disorders such asthalassemias).

4. A Complete Listing of the Equations for the Mathematical Model.

The model consists of the following structured population equationsusing the symbols defined in FIG. 9:

${{\left. {{{\left. {\left. \mspace{20mu}{{{{\frac{\partial}{\partial t}{p\left( {t,\mu^{p}} \right)}} + {\frac{\partial}{\partial\mu^{p}}{p\left( {t,\mu^{p}} \right)}}} = {\beta^{p}{p\left( {t,\mu^{p}} \right)}}},\mspace{20mu}{{{\frac{\partial}{\partial t}{q\left( {t,\mu^{q}} \right)}} + {\frac{\partial}{\partial\mu^{q}}{q\left( {t,\mu^{q}} \right)}}} = {\left( {\beta^{q} - {\alpha^{q}\left( {E(t)} \right)}} \right){q\left( {t,\mu^{q}} \right)}}},{{\frac{\partial}{\partial t}{r\left( {t,\mu^{r},h^{r},\tau^{r}} \right)}} + {\frac{\partial}{\partial\mu^{r}}{r\left( {t,\mu^{r},h^{r},\tau^{r}} \right)}} + {{v_{2}^{r}\left( {\tau,{a_{1}(t)}} \right)}\frac{\partial}{\partial h^{r}}{r\left( {t,\mu^{r},h^{r},\tau^{r}} \right)}}},{{{+ \frac{\partial}{\partial\tau^{r}}}\left( {{v_{3}^{r}\left( {\mu^{r},h^{r},\tau^{r}} \right)}{r\left( {t,\mu^{r},h^{r},\tau^{r}} \right)}} \right)} = {{8\beta\;{u\left( {t,\mu^{r},{2\; h^{r}},{2\tau^{r}}} \right)}} - {\left( {\beta + {\alpha^{r}\left( {\mu^{r},h^{r}} \right)}} \right)r\left( {t,\mu^{r},h^{r},\tau^{r}} \right)}}},{\frac{\partial}{\partial t}{s\left( {t,\mu^{s},h^{s},\tau^{s}} \right)}}} \right) + {{v_{1}^{s}\left( {E(t)} \right)}\frac{\partial}{\partial\mu^{s}}{s\left( {t,\mu^{r},h^{r},\tau^{r}} \right)}} + {{v_{2}^{s}\left( {\tau^{s},{a_{1}(t)}} \right)}\frac{\partial}{\partial h^{s}}{s\left( {t,\mu^{s},h^{s},\tau^{s}} \right)}}} \right) + {\frac{\partial}{\partial\tau^{s}}\left( {{v_{3}^{s}\left( {\mu^{s},h^{s},\tau^{s}} \right)}{s\left( {t,\mu^{s},h^{s},\tau^{s}} \right)}} \right)}} = {{- {\alpha^{3}\left( {\mu^{s},h^{s}} \right)}}{s\left( {t,\mu^{s},h^{s},\tau^{s}} \right)}}},{\frac{\partial}{\partial t}{m\left( {t,\mu^{m},h^{m},\tau^{m}} \right)}}} \right) + {\frac{\partial}{\partial\mu^{m}}{m\left( {t,\mu^{m},h^{m}} \right)}} + {{v_{2}^{m}\left( \mu_{m} \right)}\frac{\partial}{\partial h^{m}}{m\left( {t,\mu^{m},h^{m}} \right)}}} = {{- {\alpha^{m}\left( {{E(t)},\mu^{m},h^{m}} \right)}}{m\left( {t,\mu^{m}} \right)}}},\mspace{20mu}{{\frac{\mathbb{d}}{\mathbb{d}t}{E^{end}(t)}} = {{\frac{1}{TBV}{E_{in}^{end}(t)}} - {c_{\deg}^{end}{E^{end}(t)}}}},\mspace{20mu}{{\frac{\mathbb{d}}{\mathbb{d}t}{E^{ex}(t)}} = {{\frac{1}{TBV}{E_{in}^{ex}(t)}} - {c_{\deg}^{ex}{{E^{ex}(t)}.}}}}$Further, the iron model can be described as follows:

${{\frac{\mathbb{d}}{\mathbb{d}t}{a_{1}(t)}} = {{{{k_{41}\left( {H(t)} \right)}{a_{4}(t)}} + {{a_{{iv},{total}}(t)}{\delta_{t_{0}}(t)}} + {{k_{{gastro},1}\left( {H(t)} \right)}{a_{gastro}(t)}} + {{k_{51}\left( {H(t)} \right)}{a_{5}(t)}} - {k_{15}{a_{1}(t)}} - {{d_{0}\left( {{\int_{\mu_{\min}^{r}}^{\mu_{\max}^{r}}{\int_{h_{\min}^{r}}^{h_{\max}^{r}}{\int_{\tau_{\min}^{r}}^{\tau_{\max}^{r}}{\tau^{r}{r\left( {t,\mu^{r},h^{r}} \right)}\ {\mathbb{d}\mu^{r}}{\mathbb{d}h^{r}}{\mathbb{d}\tau^{r}}}}}} + {\int_{\mu_{\min}^{s}}^{\mu_{\max}^{s}}{\int_{h_{\min}^{s}}^{h_{\max}^{s}}{\int_{\tau_{\min}^{s}}^{\tau_{\max}^{s}}{\tau^{s}{s\left( {t,\mu^{s},h^{s}} \right)}\ {\mathbb{d}\mu^{s}}{\mathbb{d}h^{s}}{\mathbb{d}\tau^{s}}}}}}} \right)}{{a_{1}(t)}.{a_{2}(t)}}}} = {{{\int_{\mu_{\min}^{r}}^{\mu_{\max}^{r}}{\int_{h_{\min}^{r}}^{h_{\max}^{r}}{\int_{\tau_{\min}^{r}}^{\tau_{\max}^{r}}{h^{r}{r\left( {t,\mu^{r},h^{r}} \right)}\ {\mathbb{d}\mu^{r}}{\mathbb{d}h^{r}}{\mathbb{d}\tau^{r}}}}}} + {\int_{\mu_{\min}^{s}}^{\mu_{\max}^{s}}{\int_{h_{\min}^{s}}^{h_{\max}^{s}}{\int_{\tau_{\min}^{s}}^{\tau_{\max}^{s}}{h^{s}{s\left( {t,\mu^{s},h^{s}} \right)}\ {\mathbb{d}\mu^{s}}{\mathbb{d}h^{s}}{{\mathbb{d}\tau^{s}}.\mspace{20mu}{a_{3}(t)}}}}}}} = {{\int_{\mu_{\min}^{m}}^{\mu_{\max}^{m}}{\int_{h_{\min}^{m}}^{h_{\max}^{m}}{h^{m}{m\left( {t,\mu^{m},h^{m}} \right)}\ {\mathbb{d}h^{m}}\ {{\mathbb{d}\mu^{m}}.\frac{\mathbb{d}}{\mathbb{d}t}}{a_{4}(t)}}}} = {{\int_{\mu_{\min}^{m}}^{\mu_{\max}^{m}}{\int_{h_{\min}^{m}}^{h_{\max}^{m}}{h^{m}{\alpha^{m}\left( {{E(t)},\mu^{m},h^{m}} \right)}{m\left( {t,\mu^{m},h^{m}} \right)}\ {\mathbb{d}\mu^{m}}\ {\mathbb{d}h^{m}}}}} + {\int_{\mu_{\min}^{m}}^{\mu_{\max}^{m}}{\int_{h_{\min}^{m}}^{h_{\max}^{m}}{{v_{2}^{m}\left( \mu^{m} \right)}h^{m}{m\left( {t,\mu^{m},h^{m}} \right)}\ {\mathbb{d}\mu^{m}}\ {\mathbb{d}h^{m}}}}} + {\int_{\mu_{\min}^{r}}^{\mu_{\max}^{r}}{\int_{h_{\min}^{r}}^{h_{\max}^{r}}{\int_{\tau_{\min}^{r}}^{\tau_{\max}^{r}}{h^{r}{r\left( {t,\mu^{r},h^{r}} \right)}\ {\mathbb{d}\mu^{r}}{\mathbb{d}h^{r}}{\mathbb{d}\tau^{r}}}}}} + {\int_{\mu_{\min}^{s}}^{\mu_{\max}^{s}}{\int_{h_{\min}^{s}}^{h_{\max}^{s}}{\int_{\tau_{\min}^{s}}^{\tau_{\max}^{s}}{h^{s}{s\left( {t,\mu^{s},h^{s}} \right)}\ {\mathbb{d}\mu^{s}}{\mathbb{d}h^{s}}{\mathbb{d}\tau^{s}}}}}} - {{k_{41}\left( {H(t)} \right)}{a_{4}(t)}} - {k_{45}{a_{4}(t)}}}}}}},\mspace{20mu}{{\frac{\mathbb{d}}{\mathbb{d}t}{a_{5}(t)}} = {{k_{15}{a_{1}(t)}} + {k_{45}{a_{4}(t)}} - {{k_{51}\left( {H(t)} \right)}{a_{5}(t)}}}},\mspace{20mu}{{\frac{\mathbb{d}}{\mathbb{d}t}{H(t)}} = {{\frac{1}{TBV}{H_{prod}(t)}} - {c_{\deg}^{H}{{H(t)}.}}}}$The boundary conditions are

$\mspace{20mu}{{{p\left( {t,0} \right)} = {S_{0}\left( {E(t)} \right)}},\mspace{20mu}{{q\left( {t,\mu_{\min}^{q}} \right)} = {p\left( {t,\mu_{\max}^{p}} \right)}},\mspace{20mu}{{r\left( {t,\mu_{\min}^{r},0,\tau^{r}} \right)} = {{q\left( {t,\mu_{\max}} \right)}{\eta\left( \tau^{r} \right)}}},{{r\left( {t,\mu^{r},h^{r},\tau^{r}} \right)} = 0},{{{if}\mspace{14mu}\left( {\mu^{r},h^{r},\tau^{r}} \right){is}\mspace{14mu} a\mspace{14mu}{boundary}\mspace{14mu}{point}\mspace{14mu}{and}\mspace{14mu} h^{r}} > 0},\mspace{20mu}{{{v^{s}\left( {E(t)} \right)}{s\left( {t,\mu_{\min}^{s},h^{s},\tau^{s}} \right)}} = {r\left( {t,\mu_{\max}^{r},h^{r},\tau^{r}} \right)}},\mspace{20mu}{{m\left( {t,\mu^{m},h_{\min}^{m}} \right)} = 0},{{m\left( {t,0,h^{m}} \right)} = \left\{ \begin{matrix}{{v^{s}\left( {E(t)} \right)}{\int_{\tau_{\min}^{s}}^{\tau_{\max}^{s}}{{s\left( {t,\mu_{\max}^{s},h^{s},\tau^{s}} \right)}\ {\mathbb{d}\tau^{s}}}}} & {{{{for}\mspace{14mu}{E(t)}} \leq E^{*}},} \\{\int_{\tau_{\min}^{r}}^{\tau_{\max}^{r}}{{r\left( {t,\mu_{\max}^{r},h^{r},\tau^{r}} \right)}\ {\mathbb{d}\tau^{r}}}} & {{{{for}\mspace{14mu}{E(t)}} > E^{*}},}\end{matrix} \right.}}$and the initial values are given by

-   -   p(0,μ^(p))=p₀(μ^(p)),    -   q(0,μ^(q))=q₀(μ^(q)),    -   r(0,μ^(r),h^(r),τ^(r))=r₀(μ^(r),h^(r),τ²),    -   s(0,μ^(s),h^(s),τ^(s))=s₀(μ^(s),h^(s),τ^(s)),    -   m(0,μ^(m),h^(m))=m₀(μ^(m),h^(m)),    -   E^(end)(0)=E₀ ^(end),    -   E^(ex)(0)=E₀ ^(ex),    -   a₁(0)=a_(1,0),    -   a₂(0)=a_(2,0),    -   a₃(0)=a_(3,0),    -   a₄(0)=a_(4,0),    -   a₅(0)=a_(5,0),    -   H(0)=H₀.        Moreover, there are a number of auxiliary equations for the cell        population model

$\mspace{20mu}{{{S_{0}\left( {E(t)} \right)} = {\frac{a_{1} - b_{1}}{1 + {\mathbb{e}}^{{{- k_{1}}{E{(t)}}} + c_{1}}} + b_{1}}},\mspace{20mu}{{\alpha^{q}\left( {E(t)} \right)} = {\frac{a_{2} - b_{2}}{1 + {\mathbb{e}}^{{k_{2}{E{(t)}}} - c_{2}}} + b_{2}}},\mspace{20mu}{{v_{2}^{r}\left( {\tau^{r},a_{1}} \right)} = {d_{2}\tau^{r}a_{1}}},v_{3}^{r},{\left( {\mu^{r},h^{r},\tau^{r}} \right) = {{a_{3}{\tau^{r}\left( {h_{*}^{r} - {b_{3}\left( h^{r} \right)}^{2}} \right)}} - {c_{3}{\max\left( {\tau^{r},0} \right)}{\left( {\mu^{r} - \mu_{\max}^{q}} \right)/\left( {1 + h^{r}} \right)}}}},\mspace{20mu}{{\alpha^{r}\left( {\mu^{r},h^{r}} \right)} = {{b_{4}\left( {h^{**} - h^{r}} \right)}^{2}/\left( {1 + \left( {\mu^{**} - \left( {\mu^{r} - \mu_{\max}^{q}} \right)} \right)^{2}} \right)}},\mspace{20mu}{{\int_{\tau_{\min}}^{\tau_{\max}}{{\eta(\tau)}\ {\mathbb{d}\tau}}} = 1},\mspace{20mu}{{v_{1}^{s}\left( {E(t)} \right)} = {\frac{a_{5} - b_{5}}{1 + {\mathbb{e}}^{{{- k_{5}}{E{(t)}}} + c_{5}}} + b_{5}}},\mspace{20mu}{{v_{2}^{s}\left( {\tau_{s},a_{1}} \right)} = {d_{2}\tau^{s}a_{1}}},v_{3}^{s},{\left( {\mu^{s},h^{s},\tau^{s}} \right) = {{a_{3}{\tau^{s}\left( {h_{*}^{s} - {b_{3}\left( h^{s} \right)}^{2}} \right)}} - {c_{3}{\max\left( {\tau^{s},0} \right)}{\left( {\mu^{s} - \mu_{\max}^{r}} \right)/\left( {1 + h^{s}} \right)}}}},\mspace{20mu}{{\alpha^{s}\left( {\mu^{s},h^{s}} \right)} = {{b_{4}\left( {h^{**} - h^{s}} \right)}^{2}/\left( {1 + \left( {\mu^{**} - \left( {\mu^{s} - \mu_{\max}^{r}} \right)} \right)^{2}} \right)}},\mspace{20mu}{{v_{2}^{m}\left( \mu^{m} \right)} = {{b_{6}\left( {\mu^{m} - \mu_{*}^{m}} \right)}^{2} - c_{6}}},{{\alpha^{m}\left( {{E(t)},\mu^{m},h^{m}} \right)} = \left\{ {{\begin{matrix}{{\gamma^{m}\left( {\mu^{m},h^{m}} \right)} + {\min\left( \frac{c_{E}}{{E(t)^{k_{E}}},b_{E}} \right)}} & {{for}\mspace{14mu}\begin{matrix}{{{E(t)} < \tau_{E}},} \\{{\mu_{\min}^{m,n} \leq \mu^{m} \leq \mu_{\max}^{m,n}},}\end{matrix}} \\{\gamma^{m}\left( {\mu^{m},h^{m}} \right)} & {{otherwise},}\end{matrix}{\gamma^{m}\left( {\mu^{m},h^{m}} \right)}} = \left\{ {{{\begin{matrix}\alpha_{r}^{m} & {{for}\mspace{14mu}\begin{matrix}{{\mu^{m} \in \left\lbrack {\mu_{\min}^{m},{\mu_{\max}^{m} - \delta}} \right\rbrack},} \\{{h^{m} \in \left\lbrack {{h_{\min}^{m} + \delta},h_{\max}^{m}} \right\rbrack},}\end{matrix}} \\{{\gamma_{\mu}^{m}\left( \mu^{m} \right)}{\gamma_{h}^{m}\left( h^{m} \right)}} & {{for}\mspace{14mu}\begin{matrix}{{\mu^{m} \in \left\lbrack {{\mu_{\max}^{m} - \delta},\mu_{\max}^{m}} \right\rbrack},} \\{{h^{m} \in \left\lbrack {h_{\min}^{m},{h_{\min}^{m} + \delta}} \right\rbrack},}\end{matrix}} \\\infty & {{{{for}\mspace{14mu}\mu^{m}} \geq \mu_{\max}^{m}},{h^{m} \leq h_{\min}^{m}},}\end{matrix}\mspace{20mu}{\gamma_{\mu}^{m}\left( \mu^{m} \right)}} = \frac{3\alpha_{rand}^{m}\delta^{2}}{\left( \mu^{m} \right)^{2} - {2\left( {\mu_{\max}^{m} + \delta} \right)\mu^{m}} + {\left( {\mu_{\max}^{m} + {2\delta}} \right)\mu_{\max}^{m}}}},\mspace{20mu}{{\gamma_{h}^{m}\left( h^{m} \right)} = \frac{3\delta^{2}}{\left( h^{m} \right)^{2} - {2\left( {h_{\min}^{m} - \delta} \right)h^{m}} + {\left( {h_{\min}^{m} - {2\delta}} \right)h_{\max}^{m}}}},\mspace{20mu}{{E_{in}^{end}(t)} = {\frac{a_{7} - b_{7}}{1 + {\mathbb{e}}^{{k_{7}{{{Hb}{(t)}}/{TBV}}} - c_{7}}} + b_{7}}},\mspace{20mu}{{{Hb}(t)} = {\int_{h_{\min}^{m}}^{h_{\max}^{m}}{{m\left( {t,\mu^{m},h^{m}} \right)}\ {\mathbb{d}h^{m}}}}},\mspace{20mu}{{E(t)} = {{E^{ex}(t)} + {E^{in}(t)}}},} \right.} \right.}}$and also some auxiliary equations for the iron compartment model

$\mspace{20mu}{{{k_{{gastro},1}\left( {H(t)} \right)} = {\frac{a_{8} - b_{8}}{1 + {\mathbb{e}}^{{{- k_{g}}{H{(t)}}} + c_{8}}} + b_{8}}},\mspace{20mu}{{H_{prod}(t)} = {\frac{a_{H} - b_{H}}{1 + {\mathbb{e}}^{{{- k_{H}}{C{(t)}}} + c_{H}}} + b_{H}}},\mspace{20mu}{{C(t)} = {{{- c_{U}}{U(t)}} + {c_{a_{5}}{a_{5}(t)}} + {c_{i}{i(t)}}}},{{U(t)} = {{\int_{\mu_{\min}^{r}}^{\mu_{\max}^{r}}{\int_{h_{\min}^{r}}^{h_{\max}^{r}}{\int_{\tau_{\min}^{r}}^{\tau_{\max}^{r}}{{r\left( {t,\mu^{r},h^{r}} \right)}\ {\mathbb{d}\mu^{r}}{\mathbb{d}h^{r}}{\mathbb{d}\tau^{r}}}}}} + {\int_{\mu_{\min}^{s}}^{\mu_{\max}^{s}}{\int_{h_{\min}^{s}}^{h_{\max}^{s}}{\int_{\tau_{\min}^{s}}^{\tau_{\max}^{s}}{{s\left( {t,\mu^{s},h^{s}} \right)}\ {\mathbb{d}\mu^{s}}{\mathbb{d}h^{s}}{\mathbb{d}\tau^{s}}}}}}}},\mspace{20mu}{{k_{41}\left( {H(t)} \right)} = {\frac{a_{9} - b_{9}}{1 + {\mathbb{e}}^{{{- k_{9}}{H{(t)}}} + c_{9}}} + b_{9}}},{and}}$$\mspace{20mu}{{k_{51}\left( {H(t)} \right)} = {\frac{a_{10} - b_{10}}{1 + {\mathbb{e}}^{{{- k_{10}}{H{(t)}}} + c_{10}}} + {b_{10}.}}}$Exemplification of the Above Described Models5. Simulations for the Iron Model.

TABLE 2 Iron compartments in normal man (estimates for an average man.Height: 177 cm, weight: 70 kg. See Williams Hematology Table 40.1)Compartment Iron content (mg) Total body iron (%) Hemoglobin iron 220067 (Precursors + Erythrocytes) Storage iron 1000 27 Labile pool(macrophages) 80 2.2 Transport iron (plasma) 3 0.08 Other tissue iron138 3.72

TABLE 3 Parameter values for the iron model in a steady state. ParameterMeaning Value k₁₂ transfer rate from plasma to precursor cells 8.5 k₁₅transfer rate from plasma to storage 1.62 k₂₃ transfer rate fromprecursor cells to 0.1199 erythrocytes k₂₄ transfer rate from precursorcells to 0.0101 macrophages k₃₄ transfer rate from erythrocytes tomacrophages 0.0012 k₄₁ transfer rate from macrophages to plasma 0.29965k₄₅ transfer rate from macrophages to plasma 0.0256 k₅₁ transfer ratefrom storage to plasma 0.007

Values for the parameters of the iron submodel are assigned as describedbelow. First, parameters are formed that are able to reflect thesituation in a steady state. For this purpose, the population model andthe iron model are decoded and the compartment model describing ironhomeostasis is the only model investigated. Further, because of steadystate, all transfer rates are assumed to be constant. Then, values areassigned to the transfer rates, such that the flow of iron from onecompartment to another, and also the total amount of iron in eachcompartment, are assumed to be within physiological ranges. For adescription of the amount of normal daily iron transfer see FIG. 4.Moreover, for instance in Williams Hematology, a table can be foundwhere the iron content for a 70 kg man is presented. See Table 2 for alist of the values.

Parameter values chosen based on this information are shown in Table 3.Note that in the steady state, losses of iron and absorption of iron areequal and, therefore, are not considered here. Thus, for the moment,iron homeostasis is considered as a closed system. Using these parametervalues, a simulation was started close to the steady state values foreach compartment, as shown in FIG. 10. In addition, a situation wasconsidered where more iron is incorporated in hemoglobin and, as aconsequence, storage iron is decreased, as shown in FIG. 11. Since, inthis case the model is run with constant transfer rates and no feedbackor regulation is taken into account, it takes quite a long time untilthe system finally reaches its steady state.

Further Development of the Erythropoiesis and Iron Hemeostasis Models

The models presented in Sections 2 and 3 above are very close tophysiology and incorporate a lot of physiological mechanisms involved inerythropoiesis and iron homeostasis. As a consequence, the models arevery extensive and lead to a class of extremely complicated mathematicalequations. To solve the models, special numerical schemes have beendeveloped to find numerical approximations for the population equations.See D. H. Fuertinger. A model for erythropoiesis. PhD thesis, Universityof Graz, Austria, 2012 (hereinafter “Fuertinger Thesis”). Although theapproach chosen is very efficient, it is extremely computationallyintensive, because systems of several thousands ordinary differentialequations (ODE) have to be solved at every time step. Therefore, at thepresent time, simulations are only realizable on high performancecomputers, and thus not really practicable or viable in every dayclinical practice and patient care, where such computers are nottypically available. Since the iron component is such an important partof treatment of the anemia of renal failure, a simplified version of themodels presented in Sections 2 and 3 is described hereinafter. The modelpresented hereinafter is still comprehensive and demanding with regardto computational costs, although the complexity of the numericalapproximations reduces significantly. Instead of solving systems ofseveral thousand ODEs, less than about a hundred ODEs have to be solvedat every time step, which reduces the computational effort tremendously.

The first step is to restrict the population model presented in Section2 to population equations with only one structural attribute, namelycell age (also called maturity of the cell). Avoidance of the structuralattributes hemoglobin and transferrin receptors makes it necessary totake a step backward and reconsider how an iron model can be linked to amodel describing the erythroid cell populations. Several of theassumptions made in Sections 2 and 3 have to be revisited. It is stillimportant to keep track of how much iron is taken up by precursor cellsand how much hemoglobin is contained in every cell. The lack of thestructural feature hemoglobin, which allowed one to compute thisquantity relatively easily, requires that other mechanisms have to bedeveloped to replace it. The equations and assumptions are chosen suchthat the simplified model is as close to the one previously shown aspossible. Nevertheless, these ‘workarounds’ only enable the descriptionof the phenomenon of the accumulation of hemoglobin in cells and not somuch the physiological mechanisms involved. Whereas this might not makea difference in most cases, the more extensive model might be needed toexplain the development of RBC comprehensively for some patients withsevere changes in these mechanisms.

To ensure that the description below can be read independently, a fewimportant aspects of erythropoiesis and iron homeostasis are repeatedhereinafter.

6. Erythropoiesis Model.

The aim is to develop a mathematical model which incorporates the basicmechanisms governing erythropoiesis including the control actions oferythropoietin and the potentially restricting effect of iron. Thus, themodel has to include the regulatory processes for iron homeostasis, inparticular describing those which are connected to erythropoiesis. Theresult is a comprehensive mathematical model which is able to predictred blood cell mass under EPO- and iron-therapy.

As a first step, in this chapter a model is developed which focusessolely on the effects of erythropoietin on erythroid cells. Hence,throughout this section, erythropoiesis is assumed to be sufficientlysupplied with iron. Hereinafter, this model will be referred to as theErythropoiesis Model, described in PCT Application No.PCT/US/2012/054264 published as WO 2013/036836 A2 on Mar. 14, 2013,whose contents and teachings are incorporated herein by reference intheir entirety. The mathematical model developed in this section isapplicable for a number of different situations, as long as it isreasonable to assume that erythropoiesis is not impaired because of alack of iron availability. For instance, the Erythropoiesis Model isable to describe the recovery of red cell mass after blood donation, thereaction of the body to presurgical administration of ESAs, and changesin the number of erythrocytes of high altitude dwellers descending tosea level. See Fuertinger Thesis; D. H. Fuertinger, F. Kappel, S.Thijssen, N. W. Levin, and P. Kotanko. Journal of Mathematical Biology,66(6):1209-1240, 2013. However, the Erythropoiesis Model is onlyapplicable for a subpopulation of dialysis patients. 80-90% of dialysispatients treated with ESA suffer from functional iron deficiency at somestage in their therapy. See R. M. Schaefer and L. Schaefer. NephrologyDialysis Transplantation, 13:9-12, 1998. It has to be consideredcarefully whether, even when treated with iron compounds, it isreasonable to assume that iron supply is sufficient in patientssuffering from chronic kidney disease (CKD patients on dialysistreatment or not). Moreover, the Erythropoiesis Model helps tounderstand what the most important dynamics considering red blood cellproduction are. In Section 7 an Iron Model (compartment model) isdeveloped and in Section 7.2 it is explained how this model is linked tothe Erythropoiesis Model and extends it to a even more extensive anddetailed model for erythropoiesis, thus providing a solid basis for amore general description of the altered processes during the verycomplex pathological situation of anemia of renal failure.

The Erythropoiesis Model is based on structured population modelsdescribing the different erythroid cell stages. Five differentpopulation classes of cells are considered: BFU-Es, CFU-Es,erythroblasts, marrow reticulocytes and mature erythrocytes circulatingin the bloodstream (including peripheral reticulocytes). In theErythropoiesis Model individual cells are distinguished according totheir maturity, which can also be referred to as cell age. Thecommitment to the erythroid lineage is an irreversible event. Adifferentiated cell cannot regress or switch into anotherdifferentiation pathway. Thus, once a multipotent stem cell is committedto the erythroid lineage, it undergoes the complete series ofdifferentiations until it becomes a mature red blood cell, or it dieseventually during this process. While maturing, the cell divides anumber of times. Hence, age-structured population models of the form

$\begin{matrix}{{{\frac{\partial}{\partial t}{u\left( {t,\mu} \right)}} + {{v\left( {E(t)} \right)}\frac{\partial}{\partial\mu}{u\left( {t,\mu} \right)}}} = {\left( {{\beta\left( {{E(t)},\mu} \right)} - {\alpha\left( {{E(t)},\mu} \right)}} \right){u\left( {t,\mu} \right)}}} & (10)\end{matrix}$are used in order to describe the development of the cell populations.Here u(t,μ) denotes the population density of the cell population attime t with maturity μ. Further, the functions β and α describe thedivision rate and the rate of apoptosis of the cells, respectively. Thefunctions α and β a priori depend on the maturity μ and theconcentration of EPO E(t), respectively, at time t. The function vdescribes the maturation velocity and depends on the concentration E(t).

For the different population classes, the characteristic properties(proliferation rate, rate of apoptosis and maturation velocity) differ.While the cell matures, it changes its morphological characteristics,such as the number of EPO- and transferrin-receptors expressed on thesurface.

In the following sections, the assumptions that were made for thedifferent population classes are listed, briefly described, and themathematical equations that arise in consequence are stated. FIG. 6Agives an overview of the design and organization of the model.

6.1. Progenitor Cells: BFU-E and CFU-E Cells

Assumptions:

1. The number of cells, which commit to the erythroid lineage, isconstant.

2. Cells normally stay in this stage for 8 days (3 days BFU-E and 5 daysCFU-E).

3. EPO has no effect on the number of divisions or the rate of apoptosisof BFU-E.

4. The proliferation rate of CFU-E cells is constant.

5. The rate of apoptosis of CFU-E depends highly on EPO levels.

6. The maturation velocities of BFU-E and CFU-E cells are constant.

The process by which stem cells are recruited into proliferatingprogenitor population remains unclear. There are several hypothesesincluding an environmental dependency, that it is a random event, etc.For the moment, the number of stem cells entering the erythroid lineageis assumed to be independent of EPO and thus constant. The change inpopulation of the progenitor cells over time are described consideringtwo different classes of cells: namely BFU-E and CFU-E cells.

The earliest identifiable erythroid progenitor cell is the Burst-formingUnit Erythroid (BFU-E). At first these cells express only a very smallnumber of EPO receptors on the surface. See Wintrobe's ClinicalHematology. Thus, it is reasonable to assume EPO has no effect onproliferation or apoptosis of these cells. See Wu et al. In culture itlasts around 6-7 days until human BFU-E have all the functionalcharacteristics of the next cell stage—Colony-forming Unit Erythroidcells (CFU-E cells) (Assumption 2). See Williams Hematology.Morphologically, it is difficult to distinguish between those two typesof cells, because there are cells in between these two developmentalstages which show characteristic properties between BFU-E and CFU-Ecells. Therefore, a distinction is valid but artificial.

The BFU-E cell class is described by the following population equation

$\quad\left\{ \begin{matrix}{{\frac{\partial}{\partial t}{p\left( {t,\mu^{p}} \right)}} + {\frac{\partial}{\partial\mu^{p}}{p\left( {t,\mu^{p}} \right)}}} & {{= {\beta^{p}{p\left( {t,\mu^{p}} \right)}}},} \\{p\left( {t,0} \right)} & {{= S_{0}},} \\{p\left( {0,\mu^{p}} \right)} & {{= {p_{0}\left( \mu^{p} \right)}},}\end{matrix} \right.$

where p(t,μ^(p)) is the population density of the cell class at time twith maturity μ^(p), 0≦μ^(p)≦μ_(max) ^(p)=7, t>0. Further, β^(p) is aconstant proliferation rate and α^(p)≡0 (Assumption 3), S₀ describes thenumber of cells committing to the erythroid lineage (Assumption 1) andp₀(μ^(p)) is the population density at t=0.

Once a cell reaches the maximum age for BFU-E cells, it leaves thispopulation class and enters the CFU-E class. Consequently, there is acontinual flux of cells from one population class to the next one. CFU-Eare more rapidly dividing cells than BFU-E. See Wu et al. During thisstage, cells are very sensitive to EPO levels, and, under normalconditions, large numbers of generated CFU-E are not surviving. SeeWintrobe's Clinical Hematology. CFU-E are highly dependent on EPO toprevent them from apoptosis, i.e., the mortality for this populationclass depends on the EPO concentration. Altogether, the followingequations for the second class are obtained:

$\quad\left\{ \begin{matrix}{{\frac{\partial}{\partial t}{q\left( {t,\mu^{q}} \right)}} + {\frac{\partial}{\partial\mu^{q}}{q\left( {t,\mu^{q}} \right)}}} & {{= {\left( {\beta^{q} - {\alpha^{q}\left( {E(t)} \right)}} \right){q\left( {t,\mu^{q}} \right)}}},} \\{q\left( {t,\mu_{\min}^{q}} \right)} & {{= {p\left( {t,\mu_{\max}^{p}} \right)}},} \\{q\left( {0,\mu^{q}} \right)} & {{= {q_{0}\left( \mu^{q} \right)}},}\end{matrix} \right.$

where q(t,μ^(q)) is the population density of the CFU-E class at time twith maturity μ^(q), t>0 and 7=μ_(min) ^(q)≦μ^(q)≦μ_(max) ^(q)=13.Further on, β^(q) stands for a constant proliferation rate (Assumption4), α^(q)(E(t)) denotes the apoptosis rate depending on theEPO-concentration (Assumption 5), q(t,μ_(min) ^(q))=p(t,μ_(max) ^(p))describes the number of cells leaving the BFU-E cell stage and enteringthe CFU-E cell stage and q₀(μ^(q)) is the population density at t=0.

A sigmoid function is used to describe the rate of apoptosis,

${{\alpha^{q}\left( {E(t)} \right)} = {\frac{a_{1} - b_{1}}{1 + {\mathbb{e}}^{{k_{1}{E{(t)}}} - c_{1}}} + b_{1}}},$

where E(t) is the EPO concentration at time t and a₁, b₁, c₁ and k₁ arepositive constants with a₁>b₁. The function α^(q) monotonicallydecreases with increasing EPO concentration. Thus, a higher level of EPOcauses more cells to survive.

Note that, because of Assumption 6, the cell age of progenitor cellscoincides with the actual age of the cell, i.e., one definesv^(p)=v^(q)≡1, ∀t≧0.

6.2. Precursor Cells: Erythroblasts and Marrow Reticulocytes

Assumptions:

7. Cells stay in this stage for 6-8 days (5 days erythroblasts and 1-3days marrow reticulocytes).

8. The class erythroblasts consists of all cell stages fromproerythroblast to orthochromatophilic erythroblast.

9. EPO has no effect on the number of divisions or the rate of apoptosisof erythroblasts. The proliferation rate of erythroblasts is assumed tobe constant.

10. The maturation velocity of erythroblasts is constant.

11. Reticulocytes mature but do not proliferate.

12. The maturation velocity of reticulocytes depends on EPO.

13. A constant portion of marrow reticulocytes is phagocytosed.

After a CFU-E differentiates to a proerythroblast, it takes aboutanother 6-8 days until the cell is released from the bone marrow intothe bloodstream. See Jandl. The various stages of maturation fromproerythroblast to orthochromatophilic erythroblast (see FIG. 1) arereferred to as erythroblasts. The cells undergo several mitoticdivisions until, at the stage of orthochromatophilic erythroblast, theylose their ability to divide and enter a maturation period. Theerythroblastic pyramids appear normal, with no evidence of additionalmitotic divisions, when production increases, i.e., the proliferation oferythroblasts is assumed to be independent of EPO levels and defined tobe constant. See Williams Hematology.

Hence, the erythroblast class can be described by the following equation

$\quad\left\{ \begin{matrix}{{\frac{\partial}{\partial t}{r\left( {t,\mu^{r}} \right)}} + {\frac{\partial}{\partial\mu^{r}}{r\left( {t,\mu^{r}} \right)}}} & {{= {\beta^{r}{r\left( {t,\mu^{r}} \right)}}},} \\{r\left( {t,\mu_{\min}^{r}} \right)} & {{= {q\left( {t,\mu_{\max}^{q}} \right)}},} \\{r\left( {0,\mu^{r}} \right)} & {{= {r_{0}\left( \mu^{r} \right)}},}\end{matrix} \right.$

where r(t,μ^(r)) is the population density of the erythroblasts at timet with maturity μ^(r), t>0 and 13=μ_(min) ^(r)≦μ^(r)≦μ_(max) ^(r)=18.Further, β^(r) is a constant proliferation rate and α^(s)≡0 (Assumption9), the maturation velocity v^(r)≡1 (Assumption 10), r(t,μ_(min)^(r))=q(t,μ_(max) ^(q)) describes the number of cells leaving the CFU-Estage and entering the erythroblasts cell stage, and r₀(μ^(r)) is thepopulation density at t=0.

The differentiating process from orthochromatophilic erythroblasts tomarrow reticulocytes involves the extrusion of the cell nucleus.Reticulocytes are not capable of cell divisions (Assumption 11), i.e.,β^(s)≡0. In the Erythropoiesis Model, one does not account for animpaired erythropoiesis due to iron deficiency. Still, even when theerythroid cells in the bone marrow are sufficiently supplied with iron,not all precursor cells survive. Some of the reticulocytes die beforethey are released to the blood stream. See G. Barosi, M. Cazzola, C.Berzuini, S. Quaglini, and M. Stefanelli. British Journal ofHaematology, 61:357-370, 1985; M. Stefanelli, D. P. Bentley, I. Cavill,and H. P. Roeser. The American Journal of Physiology, 247:842-849, 1984.This is why Assumption 13 is made.

A raised EPO concentration shortens the marrow transit time of precursorcells. If the EPO level is elevated, marrow reticulocytes are releasedprematurely. This aspect is accounted for by allowing the maturationvelocity of reticulocytes to vary depending on erythropoietinconcentration. Thus, although the maximum cell age of marrowreticulocytes is fixed, the actual transit time for the cells variesbetween 1-3 days. Hence, altogether the transit time for precursor cellsis between 6-8 days (Assumption 7).

The population equation referring to marrow reticulocytes is

$\quad\left\{ \begin{matrix}{{\frac{\partial}{\partial t}{s\left( {t,\mu^{s}} \right)}} + {{v^{s}\left( {E(t)} \right)}\frac{\partial}{\partial\mu^{s}}{s\left( {t,\mu^{s}} \right)}}} & {{= {{- \alpha^{s}}{s\left( {t,\mu^{s}} \right)}}},} \\{{v^{s}\left( {E(t)} \right)}{s\left( {t,\mu_{\min}^{s}} \right)}} & {{= {r\left( {t,\mu_{\max}^{r}} \right)}},} \\{s\left( {0,\mu^{s}} \right)} & {{= {r_{0}\left( \mu^{s} \right)}},}\end{matrix} \right.$

where s(t,μ^(s)) is the reticulocytes population density at time t withmaturity μ^(s), t>0 and 18=μ_(min) ^(s)≦μ^(s)≦μ_(max) ^(s)=20. Furtheron, v^(s)(E(t)) is the maturation velocity depending on EPO (Assumption12), α^(s) denotes the rate with which reticulocytes are phagocytosed(Assumption 13), v^(s)(E(t))s(t,μ_(min) ^(s))=r(t,μ_(max) ^(r))describes the number of cells leaving the erythroblast cell stage andentering the reticulocyte cell stage, and r₀(μ^(s)) is the populationdensity at t=0.

A sigmoid function is used to describe the changes in maturationvelocity,

${{v^{s}\left( {E(t)} \right)} = {\frac{a_{2} - b_{2}}{1 + {\mathbb{e}}^{{{- k_{2}}{E{(t)}}} + c_{2}}} + b_{2}}},$

where E(t) is the EPO concentration at time t and a₂, b₂, c₂ and k₂ arepositive constants with a₂>b₂. The maturation velocity increases with arising concentration of EPO. Note that a slower maturation velocitycauses the cells to reach the maximum cell age μ_(max) ^(s) at a laterpoint, whereas a faster maturation velocity shortens the transit time,i.e., the cells reach μ_(max) ^(s) earlier.

6.3. Erythrocytes.

Assumptions:

14. Erythrocytes and blood reticulocytes are subsumed in one class.

15. Cells mature but do not proliferate.

16. There is a fixed random daily break-down of red blood cells (not tobe confused with loss of erythrocytes due to senescence).

17. A drop in EPO concentration beneath a threshold precipitatesneocytolysis.

18. Erythrocytes with age between 14-21 days are likely to be affectedby neocytolysis.

19. Cells are phagocytosed when they reach the maximum age.

20. The maximum age of erythrocytes in a healthy persons is 120 days.

The reticulocytes are released from the bone marrow into the bloodstream and within 1-2 days they mature to erythrocytes. Circulating redblood cells have no nuclei and therefore they are not able toproliferate and they are not able to repair themselves. Thus, the lifespan of these cells is limited. In healthy adults the average life spanis about 120 days, but it can significantly shorten in some pathologies.See e.g., Jandl. If not otherwise stated, 120 days is used as themaximal age of erythrocytes for the Erythropoiesis Model.

Cells can be lost for different reasons:

-   -   due to internal or external bleeding,    -   because of random daily-breakdown,    -   because neocytolysis is triggered,    -   and, because of eryptosis of senescent cells.

Altogether, the following equation is obtained for red blood cells:

$\quad\left\{ \begin{matrix}{{\frac{\partial}{\partial t}{m\left( {t,\mu^{m}} \right)}} + {\frac{\partial}{\partial\mu^{m}}{m\left( {t,\mu^{m}} \right)}}} & {{= {{- {\alpha^{m}\left( {{E(t)},\mu^{m}} \right)}}{m\left( {t,\mu^{m}} \right)}}},} \\{m\left( {t,0} \right)} & {{= {{v^{s}\left( {E(t)} \right)}{s\left( {t,\mu_{\max}^{s}} \right)}}},} \\{m\left( {0,\mu^{m}} \right)} & {{= {m_{0}\left( \mu^{m} \right)}},}\end{matrix} \right.$where m(t,μ^(m)) is the population density for the erythrocyte class attime t with maturity μ^(m), t>0 and 0=μ_(min) ^(m)≦μ^(m)≦μ_(max)^(m)=120 (Assumption 20). Moreover, m(t,0)=v^(s)(E(t))s(t,μ_(max) ^(s))describes the number of reticulocytes entering the blood stream andm₀(μ^(m)) is the population density at t=0. In the population equationα^(m)(E(t),μ^(m)) denotes a random daily break-down and neocytolysis

${\alpha^{m}\left( {{E(t)},\mu^{m}} \right)} = \left\{ \begin{matrix}{{\gamma^{m}\left( \mu^{m} \right)} + {\min\left( {{c_{E}/{E(t)}^{k_{E}}},b_{E}} \right)}} & {{{{for}\mspace{14mu}{E(t)}} < \tau_{E}},{\mu_{\min}^{m,n} \leq \mu^{m} \leq \mu_{\max}^{m,n}},} \\{\gamma^{m}\left( \mu^{m} \right)} & {{otherwise},}\end{matrix} \right.$

where b_(E), c_(E), k_(E) are positive constants, τ_(E) is the thresholdbeneath which neocytolysis is triggered (Assumption 17), and [μ_(min)^(m,n),μ_(max) ^(m,n)]=[14,21] (Assumption 18) is the age-intervalduring which cells are affected by neocytolysis, and t≧0. Additionally,it is possible to add a further term α_(bleed) ^(m)(t) to the mortalityrate in case of bleeding. Since it is assumed that a cell isphagocytosed when it reaches its maximum age (Assumption 19), themortality rate γ^(m)(μ^(m)) needs to be chosen such thatγ^(m)(μ^(m))=α_(rand) ^(m) for μ_(min) ^(m)≦μ^(m)≦μ_(max) ^(m)−δ withδ>0 sufficiently small and ∫_(μ) _(max) _(m) _(−δ) ^(μ) ^(max) ^(m)γ^(m)(μ^(m))=∞. Here α_(rand) ^(r) is a random daily break-down(Assumption 16). A possible choice for γ^(m)(μ^(m)) is

${\gamma^{m}\left( \mu^{m} \right)} = \left\{ \begin{matrix}\alpha_{rand}^{m} & {{{{for}\mspace{14mu}\mu^{m}} \in \left\lbrack {\mu_{\min}^{m},{\mu_{\max}^{m} - \delta}} \right\rbrack},} \\\frac{3\alpha_{rand}^{m}\delta^{2}}{\left( \mu^{m} \right)^{2} - {2\left( {\mu_{\max}^{m} + \delta} \right)\mu^{m}} + {\left( {\mu_{\max}^{m} + {2\delta}} \right)\mu_{\max}^{m}}} & {{{{for}\mspace{14mu}\mu^{m}} \in \left\lbrack {{\mu_{\max}^{m} - \delta},\mu_{\max}^{m}} \right)},} \\\infty & {{{for}\mspace{14mu}\mu^{m}} \geq {\mu_{\max}^{m}.}}\end{matrix} \right.$6.4. ErythropoietinAssumptions:

21. Release of EPO is controlled by a negative feedback mechanismaccording to the oxygen content.

22. Oxygen carrying capacity is directly proportional to the number oferythrocytes.

23. The degradation rate of EPO is constant.

24. There is a slight delay in reaction of the EPO production rate tothe number of RBC but this is negligible compared to the duration ofdevelopment of erythrocytes.

The kidneys adjust the release of EPO according to the oxygen carryingcapacity of the blood. If blood oxygen content is lower than normal,then EPO production increases and vice versa. Thus, the production ofEPO is controlled by a negative feedback mechanism and allows for morered blood cells to be developed in case of an undersupply of the bodywith oxygen by opposing programmed cell death of erythroid progenitorcells. Further, if EPO decreases beneath a certain threshold,neocytolysis is triggered, that is, a process wherein macrophages startto phagocytose young erythrocytes (neocytes).

A sigmoid function which depends on blood oxygen partial pressure isused to model the feedback involving the release of erythropoietinE_(in) ^(end)(t) from the kidneys into plasma. As a consequence ofAssumption 22, the amount of E_(in) ^(end)(t) of EPO released by thekidney per unit time can be directly computed by use of the totalpopulation of erythrocytes M(t)=∫₀ ^(μ) ^(max) ^(m) m(t,μ^(m))dμ^(m).Recall that this class consists of all circulating red blood cells(Assumption 14).

${{E_{in}^{end}(t)} = {\frac{a_{3} - b_{3}}{1 + {\mathbb{e}}^{{k_{3}{\overset{\sim}{M}{(t)}}} - c_{3}}} + b_{3}}},$

where {tilde over (M)}(t)=10⁻⁸ M(t)/TBV is a scaled erythrocytes“concentration”. Here, TBV is the total blood volume. The constants a₃,b₃, c₃, k₃ are positive and satisfy a₃>b₃. The function E_(in) ^(end)(t)is monotonically decreasing. Thus, the release of EPO increases if thenumber of circulating red blood cells decreases (Assumption 21). Thedynamics of the endogenous EPO concentration E^(end)(t) in plasma aredescribed by the following ordinary differential equation:

$\begin{matrix}{{{\frac{\mathbb{d}}{\mathbb{d}t}{E^{end}(t)}} = {{\frac{1}{T\; B\; V}{E_{in}^{end}(t)}} - {c_{\deg}^{end}{E^{end}(t)}}}},} & (11)\end{matrix}$where E^(end)(t) is the endogenous EPO concentration in plasma, E_(in)^(end) is the amount of EPO released by the kidneys and c_(deg) ^(end)describes the constant degradation rate of endogenous EPO (Assumption23).

The degradation rate c_(deg) ^(ex) for exogenous EPO differs from theone for endogenous EPO and varies according to the kind of ESAadministered. Therefore, an additional ODE is needed to describe thechange in the plasma concentration of an ESA

$\begin{matrix}{{\frac{\mathbb{d}}{\mathbb{d}t}{E^{ex}(t)}} = {{\frac{1}{T\; B\; V}{E_{in}^{ex}(t)}} - {c_{\deg}^{ex}{{E^{ex}(t)}.}}}} & (12)\end{matrix}$Here E_(in) ^(ex)(t) is the rate at which the artificial hormone isadministered and c_(deg) ^(ex) is the rate with which the exogenoushormone is degraded. In intravenous administration, the total amount ofthe agent is injected into a vein, within a very short time interval. Inthis case E_(in) ^(ex)(t) can be approximated by E₀ ^(ex)(t)δ_(t) ₀ (t),where E₀ ^(ex) is the amount of artificial hormone administered andδ_(t) ₀ (t) is the Dirac delta impulse located at t₀, the time when theadministration takes place. The overall concentration of EPO in bloodconsists of the naturally produced erythropoietin in the body and theadministered ESAE(t)=E ^(ex)(t)+E ^(end)(t).6.5. Mathematical Model—Erythropoiesis Model

For the convenience of the reader, all equations of the ErythropoiesisModel are collected in this section:

$\mspace{20mu}{{{{\frac{\partial}{\partial t}{p\left( {t,\mu^{p}} \right)}} + {\frac{\partial}{\partial\mu^{p}}{p\left( {t,\mu^{p}} \right)}}} = {\beta^{p}{p\left( {t,\mu^{p}} \right)}}},\mspace{20mu}{{{\frac{\partial}{\partial t}{q\left( {t,\mu^{q}} \right)}} + {\frac{\partial}{\partial\mu^{q}}{q\left( {t,\mu^{q}} \right)}}} = {\left( {\beta^{q} - {\alpha^{q}\left( {E(t)} \right)}} \right){q\left( {t,\mu^{q}} \right)}}},\mspace{20mu}{{{\frac{\partial}{\partial t}{r\left( {t,\mu^{r}} \right)}} + {\frac{\partial}{\partial\mu^{r}}{r\left( {t,\mu^{r}} \right)}}} = {\beta^{r}{r\left( {t,\mu^{r}} \right)}}},\mspace{20mu}{{{\frac{\partial}{\partial t}{s\left( {t,\mu^{s}} \right)}} + {{v^{s}\left( {E(t)} \right)}\frac{\partial}{\partial\mu^{s}}{s\left( {t,\mu^{s}} \right)}}} = {{- \alpha^{s}}{s\left( {t,\mu^{s}} \right)}}},\mspace{20mu}{{{\frac{\partial}{\partial t}{m\left( {t,\mu^{m}} \right)}} + {\frac{\partial}{\partial\mu^{m}}{m\left( {t,\mu^{m}} \right)}}} = {{- {\alpha^{m}\left( {{E(t)},\mu^{m}} \right)}}{m\left( {t,\mu^{m}} \right)}}},\mspace{20mu}{{\frac{\mathbb{d}}{\mathbb{d}t}{E^{end}(t)}} = {{\frac{1}{T\; B\; V}{E_{in}^{end}(t)}} - {c_{\deg}^{end}{E^{end}(t)}}}},\mspace{20mu}{{\frac{\mathbb{d}}{\mathbb{d}t}{E^{ex}(t)}} = {{\frac{1}{T\; B\; V}{E_{in}^{ex}(t)}} - {c_{\deg}^{ex}{E^{ex}(t)}}}},\mspace{20mu}{where}}$$\mspace{20mu}{{{\alpha^{q}\left( {E(t)} \right)} = {\frac{a_{1} - b_{1}}{1 + {\mathbb{e}}^{{k_{1}{E{(t)}}} - c_{1}}} + b_{1}}},\mspace{20mu}{{v^{s}\left( {E(t)} \right)} = {\frac{a_{2} - b_{2}}{1 + {\mathbb{e}}^{{{- k_{2}}{E{(t)}}} - c_{2}}} + {b_{2}.\mspace{20mu}{and}}}}}$${\alpha^{m}\left( {{E(t)},\mu^{m}} \right)} = \left\{ \begin{matrix}{{{\gamma^{m}\left( \mu^{m} \right)} + {\min\left( {\frac{c_{E}}{{E(t)}^{k_{E}}},b_{E}} \right)}},} & {{{{for}\mspace{14mu}{E(t)}} < \tau_{E}},{\mu_{\min}^{m,n} \leq \mu^{m} \leq \mu_{\max}^{m,n}},} \\{{\gamma^{m}\left( \mu^{m} \right)},} & {{otherwise},}\end{matrix} \right.$with

${\gamma^{m}\left( \mu^{m} \right)} = \left\{ \begin{matrix}\alpha_{rand}^{m} & {{{{for}\mspace{14mu}\mu^{m}} \in \left\lbrack {\mu_{\min}^{m},{\mu_{\max}^{m} - \delta}} \right\rbrack},} \\\frac{3\alpha_{rand}^{m}\delta^{2}}{\left( \mu^{m} \right)^{2} - {2\left( {\mu_{\max}^{m} + \delta} \right)\mu^{m}} + {\left( {\mu_{\max}^{m} + {2\delta}} \right)\mu_{\max}^{m}}} & {{{{for}\mspace{14mu}\mu^{m}} \in \left\lbrack {{\mu_{\max}^{m} - \delta},\mu_{\max}^{m}} \right)},} \\\infty & {{{for}\mspace{14mu}\mu^{m}} \geq {\mu_{\max}^{m}.}}\end{matrix} \right.$

Further, the boundary conditions are given by

-   -   p(t,0)=S₀,    -   q(t,μ_(min) ^(q))=p(t,μ_(max) ^(p)),    -   r(t,μ_(min) ^(r))=q(t,μ_(max) ^(q)),    -   v^(s)(E(t))s(t,μ_(min) ^(s))=r(t,μ_(max) ^(r)),    -   m(t,0)=v^(s)(E(t))s(t,μ_(max) ^(s)),    -   m(0,μ^(m))=m₀(μ^(m)),

and the initial values are

-   -   p(0,μ^(p))=p₀(μ^(p)),    -   q(0,μ^(q))=q₀(μ^(q)),    -   r(0,μ^(r))=r₀(μ^(r)),    -   s(0,μ^(s))=r₀(μ^(s)),    -   E^(end)(0)=E₀ ^(end),    -   E^(ex)(0)=E₀ ^(ex).

The feedback is described by

${{E_{in}^{end}(t)} = {\frac{a_{3} - b_{3}}{1 + {\mathbb{e}}^{{k_{3}{\overset{\sim}{M}{(t)}}} - c_{3}}} + b_{3}}},$

where {tilde over (M)}(t)=10⁻⁸M(t)/TBV, with M(t)=∫₀ ^(μ) ^(max) ^(m)m(t,μ^(m))dμ^(m) and one hasE(t)=E ^(ex)(t)+E ^(end)(t).7. Iron Homeostasis.

As described above, much of the iron in the human body can be found incirculating erythrocytes incorporated in hemoglobin. However, most ofthe iron used for the production of new erythrocytes comes fromhemoglobin recycling and not from absorption. When a senescent red cellis phagocytosed, the macrophage internalizes hemoglobin, degrades it andexports it back out into blood circulation.

Iron is stored in form of ferritin or hemosiderin, where the latter ismore a kind of a long-term storage, because iron within deposits ofhemosiderin is very poorly available. Storage sites can be found mainlyin the liver, the spleen and the bone marrow but, for instance, smallamounts of hemosiderin can be found in almost any tissue. There aredifferences between men and women in the amount of iron that is stored(Females about 200-400 mg, Males about 1000 mg). See Crepaldi.

Iron is carried around in plasma by a transport protein, which is calledtransferrin. Transferrin is produced by the liver. Cells that are inneed for iron express transferrin receptors, and after transferrin hasformed a complex with this receptor, iron is transported into the cell.Transferrin exists in three different forms: iron-free (apotransferrin),one iron (monoferric transferrin) and two iron (differic transferrin)atoms bound. The concentration of iron and the amount of transferrinpresent in plasma determine the relative abundance of each form.Although iron bound to transferrin amounts to only about 3 mg of totalbody iron, it is the most important iron pool with the highest rate ofturnover, about 30 mg/day. Approximately 80% of this iron is transportedto the bone marrow and used by erythroid cells. See Lieu et al. Theamount of transferrin-bound iron is determined by three processes:macrophage iron recycling, (hepatic) iron storage and duodenal ironabsorption. See Wrighting.

In the last 15 years, understanding of the regulation of systemic ironhomeostasis has changed substantially. The antimicrobial peptidehepcidin, which is secreted by the liver, was identified to be thesystemic iron regulator. Hepcidin inhibits iron transport by binding tothe iron exporter ferroportin. Hence, it keeps gut enterocytes fromsecreting iron into circulation, thereby functionally reducing ironabsorption. Further, it prevents iron release from macrophages andblocks the stores as well, by shutting off the means of transport out ofcells. See FIG. 4 for a diagram of the iron cycle and daily turnoverrates and FIG. 5 for a sketch of systemic iron regulation.

There exist some sensing mechanisms in the liver that gauge the amountof iron needed to support erythropoiesis. Note that hepatocytes in theliver, which secrete hepcidin, also release transferrin into circulationand are one of the major iron storage sites. If the concentration ofiron stores are too low and/or erythropoiesis is increased, then thelevel of the hormone hepcidin is decreased. Thus, for instance, lowhepcidin concentrations are observed in patients with absoluteiron-deficiency anemia, or anemias with high erythropoietic activity.See Goodnough 2010. The consequence of a lower hepcidin level is thatmore iron is taken up via the duodenum and more iron is released frommacrophages and from the stores. See Crichton; Fleming. If iron storesare full and/or erythropoiesis is decreased, then more hepcidin issecreted by hepatocytes. Dietary iron absorption, macrophage ironrecycling and release of iron from stores is partly or fully blocked asa result of the binding to the iron exporter ferroportin.

Hepcidin is an acute phase-reactant, and, during inflammation, cytokinesstimulate overproduction of hepcidin. Hence, the feedback is alteredduring inflammation, which is a prominent problem in chronic kidneydisease. As a consequence, the body is not able to sufficiently supplydeveloping red blood cells with iron and erythropoiesis is impaired,although there might be more than enough iron in stores and macrophages,but it is simply not available. This paradoxical situation is calledfunctional iron deficiency, a condition where stored iron is sufficient,but circulating iron is deficient, whereas a situation where the ironcontent of the body is too low (depleted stores) is referred to asabsolute iron deficiency that occurs when iron stores are depleted.

7.1 Iron Model

Assumptions

25. Five iron pools are considered within the body: a plasma pool, ironin precursor cells, iron in erythrocytes, iron in the macrophages and astorage pool.

26. Other iron in the body is ignored.

27. Hepcidin regulates iron homeostasis.

28. Stressed erythropoiesis and depleted stores decrease the amount ofhepcidin released into plasma.

29. Large iron stores and inflammation increase hepcidin levels.

30. The amount of iron in plasma is proportional to the amount oftransferrin molecules carrying iron. (Free iron is neglected.)

31. There is no differentiation between monoferric and differictransferrin.

32. Iron is only lost via loss of cells, i.e. iron lost via urine andsweat is neglected.

Iron is a part of hemoglobin which provides the means of O₂ transport tothe tissue. The rate at which iron can be released to plasma fromexisting stores and macrophages may easily limit the rate of irondelivery for hemoglobin synthesis. Therefore, a strong relationshipbetween impaired erythropoiesis and total body iron burden exists. SeeCrichton for an excellent description of iron homeostasis and itseffects.

Iron is a substance which is very tightly controlled from the body asiron overload is toxic. Iron homeostasis can be only achieved by thecontrol of absorption because the human body, unlike other mammals, isnot able to influence the excretion of iron. Under normal circumstancesthe average daily loss of iron is very small (men 1 mg, women 2 mg,because of menstruation).

For the Iron Model, five dynamic iron pools are considered within thebody: a plasma pool, iron in precursor cells, iron in erythrocytes, ironin the macrophages and a storage pool. Hence, ferritin and hemosiderinstores are lumped together. Other iron in the body is ignored. For asketch of the organizational structure of the Iron Model see FIG. 8. Theamount of iron lost via urine and sweat is neglected. Thus, in thismodel, excretion of iron only takes place when cells are lost. On onehand, this can be due to loss of RBCs, i.e., due to external bleeding,and, on the other hand, there is a daily loss of epithelial cells.Further, one assumes that only the precursor cells in the bone marrowuptake iron, and ineffective erythropoiesis depends on the amount ofiron which can be provided for erythropoiesis. Moreover, the modelaccounts for the loss of iron from erythrocytes during senescence.

The general form of the compartmental Iron Model is as follows:

${{\frac{\mathbb{d}}{\mathbb{d}t}a_{1}} = {{{k_{41}(H)}a_{4}} + {k_{{iv},1}a_{iv}} + {{k_{{gastro},1}(H)}a_{gastro}} + {{k_{51}(H)}a_{5}} - {k_{15}a_{1}} - {k_{12}a_{1}}}},\mspace{20mu}{{\frac{\mathbb{d}}{\mathbb{d}t}a_{2}} = {{k_{12}a_{4}} - f_{24} - f_{23}}},\mspace{20mu}{{\frac{\mathbb{d}}{\mathbb{d}t}a_{3}} = {f_{23} - f_{34} - f_{3,{out}}}},\mspace{20mu}{{\frac{\mathbb{d}}{\mathbb{d}t}a_{4}} = {f_{24} + {f_{34}a_{3}} - {{k_{41}(H)}a_{4}} - {k_{45}a_{4}}}},\mspace{20mu}{{\frac{\mathbb{d}}{\mathbb{d}t}a_{5}} = {{k_{15}a_{1}} + {k_{45}a_{4}} - {{k_{51}(H)}{a_{5}.}}}}$

Note, in the above equations, the dependence of particular components ont was not explicitly stated to simplify notation. Furthermore, a_(i),i=1, . . . , 5, denotes the amount of iron in the compartment i andk_(ij), i, j=1, . . . , 5 are the transfer rates from compartment i tocompartment j, i.e., k_(ij)a_(i) is the rate at which iron is movingfrom compartment i to compartment j. Further, ƒ_(ij) denote the rates atwhich iron is transferred from compartment i into compartment j, when itis not of the simple form k_(ij)a_(i); ƒ_(3,out) is the amount of ironwhich is lost by bleeding. Additionally, a_(iv) denotes the amount ofiron intravenously administered, whereas a_(gastro) is the amount ofiron in the duodenum. Finally, k_(iv,1) and k_(gastro,1)(H) are thecorresponding flow rates. The rate k_(gastro,1) can be described by adecreasing sigmoidal function

${{k_{{gastro},1}\left( {H(t)} \right)} = {\frac{a_{5} - b_{5}}{1 + {\mathbb{e}}^{{{- k_{5}}{H{(t)}}} + c_{5}}} + b_{5}}},$

where a₅, b₅, c₅ and k₅ are positive constants with a₈>b₈. The functionH(t) is the solution of the differential equation (13) below.

Since the time needed to administer iron is very short compared to othertime constants in the system, the following can be takenk _(iv,1) a _(iv) =a _(iv,total)δ_(t) ₀ .

Here a_(iv,total) is the total amount of iron administered and δ_(t) ₀is the Dirac delta function located at t₀.

7.1.1 Hepcidin Feedback

The regulator of iron homeostasis is hepcidin, a peptide hormoneproduced by the liver. Hepcidin negatively regulates the availability ofiron in plasma by direct inhibition of ferroportin. Ferroportin is theiron exporter required for iron egress from iron exporting cells, as forinstance, enterocytes and macrophages. Hepcidin acts as a medium towithhold iron from certain invasive bacteria, as they are not able toproliferate under conditions of insufficient iron supply. The dynamicsof the hepcidin concentration in plasma is calculated according to

$\begin{matrix}{{{\frac{\mathbb{d}}{\mathbb{d}t}{H(t)}} = {{\frac{1}{T\; B\; V}{H_{prod}(t)}} - {c_{\deg}^{H}{H(t)}}}},} & (13)\end{matrix}$where H_(prod)(t) is the rate at which hepcidin, produced by thehepatocytes, is released into the plasma at time t. TBV is the totalblood volume, and c_(deg) ^(H) is the rate at which hepcidin isdegraded. The production rate H_(prod)(t) of hepcidin is obtained as theresult of a feedback loop involving the total precursor cell population,the amount of iron stored and the state of inflammation of the patient,H _(prod)(t)=ƒ(U(t),a ₅(t),i(t)).

Here

U(t) = ∫_(μ_(min)^(r))^(μ_(max)^(r))r(t, μ^(r)) 𝕕μ^(r) + ∫_(μ_(min)^(s))^(μ_(max)^(s))s(t, μ^(s)) 𝕕μ^(s),

where U(t) describes the total precursor cell population(erythroblasts+reticulocytes). Further, a₅(t) is the amount of ironstored and i(t) is the inflammation status at time t. The feedback ismodeled using a monotonically increasing sigmoidal function

${{H_{prod}(t)} = {\frac{a_{H} - b_{H}}{1 + {\mathbb{e}}^{{{- k_{H}}{C{(t)}}} + c_{H}}} + b_{H}}},$

where a_(H), b_(H), c_(H) and k_(H) are positive constants witha_(H)>b_(H). The functionC(t)=−c _(U) U(t)+c _(a) ₅ a ₅(t)+c _(i) i(t)

is an auxiliary equation that increases if stores or inflammation statusrises, and decreases when more RBCs are produced. Here, c_(U), c_(a) ₅and c_(i) are positive constants.

7.2. Iron Model Linked to the Erythropoiesis Mode. 1

Assumptions

33. All types of precursor cells synthesize hemoglobin.

34. The rate at which the amount of hemoglobin in precursor cellschanges, depends on the number of transferrin receptors (TfR) on thecell membrane and the amount of iron in plasma.

35. Cells with the same age express the same amount of TfR on themembrane.

36. An erythroblast starts with 300,000 TfR on the cell membrane.

37. This number increases sharply and peaks after 1.5 days at 800,000TfR.

38. After that time, it declines to 100,000 TfR on bone marrowreticulocytes.

39. The number of TfR on bone marrow reticulocytes is constant (100,000TfR/cell).

40. The rate of apoptosis of precursor cells (erythroblasts ANDreticulocytes) depend on cell age and hemoglobin.

41. Iron that enters a precursor cell is incorporated into hemoglobin,i.e., no free iron in the cell is considered.

42. Erythrocytes don't lose iron during senescence.

Some thought has to be given to how the iron model can be linked to thepopulation model. Iron is consumed by precursor cells (erythroblasts andbone marrow reticulocytes) and used for hemoglobin synthesis. The amountof iron that is taken up by a cell depends on the number of transferrinreceptors (TfR) expressed on the membrane and the concentration of ironin plasma. In order to simplify notation, the minimal cell age ofprecursor cells μ_(min) ^(r) is set to zero, i.e., the cell age forerythroblasts is shifted from 8=μ_(min) ^(r)≦μ^(r)≦μ_(max) ^(r)=13 to0=μ_(min) ^(r)≦μ^(r)≦μ_(max) ^(r)=5 and the cell age for bone marrowreticulocytes consequently is shifted from 13=μ_(min) ^(s)≦μ^(s)≦μ_(max)^(s)=15 to 5=μ_(min) ^(s)≦μ^(s)≦μ_(max) ^(s)=7.

The contact rate between the TfRs of a cell and transferrin molecules inplasma carrying iron is governed by the law of mass action and thus isgiven by:

$\begin{matrix}{{d_{1}{\tau(t)}\frac{a_{1}(t)}{T\; B\; V}},} & (14)\end{matrix}$where d₁ is a rate constant that differs for different patients, as itdepends on the volume of the active bone marrow. Further, τ(t) is thecumulative amount of TfR on cell membranes of precursor cells(erythroblasts and bone marrow reticulocytes) at time t and a₁(t) is theconcentration of iron in plasma at time t. Note that this is at the sametime the rate with which iron leaves the plasma compartment and entersthe precursor cells. The cumulative sum of TfR expressed on allprecursor cells τ(t) is determined by

$\begin{matrix}{{{\tau(t)} = {{\int_{\mu_{\min}^{r}}^{\mu_{\max}^{r}}{{\varphi\left( \mu^{r} \right)}\ {r\left( {t,\mu^{r}} \right)}{\mathbb{d}\mu^{r}}}} + {\int_{\mu_{\min}^{s}}^{\mu_{\max}^{s}}{{\varphi\left( \mu^{s} \right)}\ {s\left( {t,\mu^{s}} \right)}{\mathbb{d}\mu^{s}}}}}},} & (15)\end{matrix}$where φ(μ) is the function describing the number of TfR on a cell withcell age μ, r(t,μ^(r)) is the density of erythroblasts, and s(t,μ^(s))is the density of bone marrow reticulocytes.7.2.1. Function for the Transferrin Receptors.

Assumptions 35-39 govern the development of the transferrin receptors.One ensures that Assumptions 35-38 are met by choosing a polynomial ofdegree 3φ(μ)=a ₄μ³ +b ₄μ² +c ₄ μ+d ₄,

and choosing the coefficients a₄, b₄, c₄ and d₄ such that the function φsatisfies the following conditions:φ(0)=30,φ(1.5)=80,φ′(1.5)=0,φ(5)=10.  (16)For computational reasons, the number of TfR is multiplied by 10⁻⁴. Thenthe first condition in eq. (16) coincides with Assumption 36. The secondand third condition reflect Assumption 37 and, with the last condition,one ensures that there is a smooth transition in the number of TfR fromerythroblastic cells to reticulocytes. Deriving and solving the linearsystem that arises from eq. (16) gives a₄=3.3016 b₄=−32.127 d₄=30.Further, using Assumption 38, the function φ can be extended to bonemarrow reticulocytes. Altogether, one has

$\begin{matrix}{{\varphi(\mu)} = \left\{ \begin{matrix}{{{3.3016\mspace{14mu}\mu^{3}} - {32.127\mspace{14mu}\mu^{2}} + {74.0952\mspace{14mu}\mu} + 30},} & {{{{for}\mspace{14mu}\mu} \in \left\lbrack {\mu_{\min}^{r},\mu_{\max}^{r}} \right\rbrack},} \\{10,} & {{{for}\mspace{14mu}\mu} \in {\left\lbrack {\mu_{\min}^{s},\mu_{\max}^{s}} \right\rbrack.}}\end{matrix} \right.} & (17)\end{matrix}$For a plot of the function φ see FIG. 12.7.2.2. Accumulation of Hemoglobin

Red blood cells carry hundreds of hemoglobin molecules which enablesthem to carry oxygen. Hemoglobin has four subunits each containing aheme group. An iron atom is embedded in the middle of the heme group.Hence, in one hemoglobin molecule four iron atoms are incorporated. Inorder to keep track of the accumulation of hemoglobin molecules in theprecursor cells, the development of the TfR and the amount of ironavailable at certain times have to be considered. A consequence ofAssumption 35 is that cells with the same cell age all carry the sameamount of hemoglobin. The amount of hemoglobin a cell with age μaccumulated can vary over time t but is invariant among cells with age μat a specific time t. First, one focuses on the erythroblast populationclass, where one has the special case that the maturation velocity v≡1,i.e. the cell age and the time for how long a cell has been anerythroblast, and thus had time to accumulate hemoglobin, coincide.Hereinafter, the actual age of the cell, i.e. the time that elapsedsince the cell became a precursor cell, will be called calendric age.This is to emphasize the difference between the cell age and the actualage, i.e. the calendric age.

The rate with which hemoglobin increases in a cell is given by

$\begin{matrix}\begin{matrix}{\frac{\mathbb{d}}{\mathbb{d}t}\left( {{h\left( \left( {t,{\mu(t)}} \right) \right)} = {{\frac{\partial}{\partial t}{h\left( {t,{\mu(t)}} \right)}} + {\frac{\partial}{\partial\mu}{h\left( {t,{\mu(t)}} \right)}{v\left( {E(t)} \right)}}}} \right.} \\{{= {{\mathbb{d}\;{\varphi\left( {\mu(t)} \right)}}\frac{a_{1}(t)}{T\; B\; V}}},}\end{matrix} & (18)\end{matrix}$where h(t,μ(t)) denotes the amount of hemoglobin in a cell with age μ(t)at time t. Note that, in most cases, the dependence of the cell age μ ontime t is not explicitly stated, because there is no necessity to do so.Here, one has to consider this fact, as it is important for thederivation of the function h. Further, d is a rate constant which isdifferent from d₁ because here one considers hemoglobin instead of iron,see eq. (14). The ratio between d and d₁ is given by the fact that eachhemoglobin molecule consists of four iron atoms. The function φ(μ) isdetermined by eq. (17), a₁(t) describes the amount of iron in the plasmacompartment and TBV denotes the total blood volume. For a detaileddescription of how the amount of hemoglobin in a cell with age μ at timet can be calculated from eq. (18), see Appendix A4.Maturation velocity v≡1

Cells in the erythroblast population class in the Erythropoiesis Modelsatisfy the condition that cell age and calendric age of the cellcoincide, see Section 6. This makes it easier to keep track of theaccumulation of hemoglobin. In order to calculate the amount ofhemoglobin a cell with cell age μ at time t carries, one has to considerthat the cell entered the population class at time t−μ+μ_(min), whereμ_(min) denotes the minimal cell age of the population class. Further,the amount of TfR on the cell membrane of the cell that entered thepopulation at time t−μ+μ_(min) was φ(μ_(min)) and the amount of iron inplasma at that time was a₁(t−μ+μ_(min)). With these considerations, theamount of hemoglobin a cell with cell age μ at time t carries is givenby the following integral

${h\left( {\overset{\_}{t},\overset{\_}{\mu}} \right)} = {\frac{\mathbb{d}}{T\; B\; V}{\int_{\overset{\_}{t} - \overset{\_}{\mu} + \mu_{\min}}^{\overset{\_}{t}}{{\varphi\left( {s - \overset{\_}{t} + \overset{\_}{\mu}} \right)}{a_{1}(s)}\ {{\mathbb{d}s}.}}}}$

Note that the calculation makes use of the history of the availableamount of iron in plasma from t−μ+μ_(min) to t.

For the erythroblast population (the cell age was transformed such thatμ_(min) ^(r)≡0) the above equation simplifies to

$\begin{matrix}{{h\left( {t,\mu^{r}} \right)} = {\frac{\mathbb{d}}{T\; B\; V}{\int_{t - \mu^{r}}^{t}{{\varphi\left( {s - t + \mu^{r}} \right)}{a_{1}(s)}\ {{\mathbb{d}s}.}}}}} & (19)\end{matrix}$Maturation velocity v=v(E(t))

In the case that cells have a variable maturation velocity, it gets moredifficult to keep track of the accumulation of hemoglobin. First, onehas to calculate at which time {circumflex over (t)} a cell that hascell age μ at time t entered the population class. Therefore, thefollowing equation has to be solved with regard to {circumflex over (t)}

${\overset{\_}{\mu} = {\mu_{\min} + {\int_{\hat{t}}^{\overset{\_}{t}}{{v\left( {E(s)}\  \right)}{\mathbb{d}s}}}}},$

i.e., one has to (numerically) find the root of the following function

$\begin{matrix}{{f\left( {\overset{\_}{t},\overset{\_}{\mu}} \right)} = {\mu_{\min} - \overset{\_}{\mu} + {\int_{\hat{t}}^{\overset{\_}{t}}{{v\left( {E(s)} \right)}{{\mathbb{d}s}.}}}}} & (20)\end{matrix}$The value η=t−{circumflex over (t)} reflects the calendric age of thecell. (It represents the time that passed since the cell entered thepopulation class and does not necessarily coincide with the cellage/maturity of the cell.)

In the Erythropoiesis model presented in Section 6, only thereticulocyte population has a variable maturation velocity. Assumption39 states that cells in this specific population class express aconstant number of TfR on their membranes. For the sake of simplicity,the special case φ≡const. is considered here, and the general case whenφ=φ(μ) is not considered.

Thus, the amount of hemoglobin in a reticulocyte with calendric ageη^(s)ε[5,8] at time t is given

$\begin{matrix}{{h\left( {t,\eta^{s}} \right)} = {{h\left( {{t - \eta^{s} + \mu_{\min}^{s}},\mu_{\max}^{r}} \right)} + {\varphi\frac{\mathbb{d}}{T\; B\; V}{\int_{t - \eta^{s} + \mu_{\min}^{s}}^{t}{{a_{1}(s)}\ {{\mathbb{d}s}.}}}}}} & (21)\end{matrix}$The calendric age of the cell has been derived using the solution of eq.(20) and h(t−η^(s)+μ_(min) ^(s),μ_(max) ^(r)) is derived from eq. (19).7.2.3 Rate of Apoptosis of Precursor Cells

The construction of the rate of apoptosis of precursor cells is inspiredby the fact that cells that do not synthesize hemoglobin sufficientlyare more likely to die. See Schmidt.

Hereinafter, one is going to refer to cell age of the erythroblasts andthe calendric age (see eq. 20) of the reticulocytes, respectively, as η.Consequently, ηε[0,8] and denotes the time that elapsed since a cellbecame a precursor cell.

First, a “reference hemoglobin accumulation curve” denoted by h*(η) iscreated using eqs. (19) and (21). The construction is based on anaverage male healthy adult, assuming that:

-   -   the person's erythropoiesis and iron homeostasis is in steady        state,    -   the total blood volume is 5 liters,    -   and the amount of iron in plasma is 3 mg, i.e., a₁≡3.

Further, healthy persons produce normocytic red blood cells and theamount of Hgb per cell lies within 26-32 pg. The contact rate d ischosen such that, under the above assumptions and by applying eqs. (19)and (21), respectively, the hemoglobin content of a cell lies wellwithin this range when the cell leaves the reticulocyte population andenters the blood stream. FIG. 13 and FIGS. 14A-14B depict theaccumulation of hemoglobin in precursor cells for d=0.215. Note that thelabel ‘days since cell became a precursor cell’ is similar to thecalendric age η and for erythroblast (day 0-5) coincides with the cellage, whereas for the reticulocytes the cell age and the calendric agemight differ. Reticulocytes may stay in the bone marrow 1-3 days, butthey all leave the bone marrow when they reach maximal maturity, i.e.,when their cell age turns 2. This means that reticulocytes might havedifferent timespans to gather iron and synthesize hemoglobin—dependingon the maturation velocity v(E(t)), see Section 6.

Next, the reference curve is shifted, by slightly reducing orincreasing, respectively, the amount of iron in plasma (a₁≡2.75 mg anda₁≡3.25 mg, respectively). The consequence being that a corridordevelops around the reference curve. In FIG. 15 the reference curve isshown as a solid line and the shifted curves appear as dashed lines. Thenext step is to design functions in a way that cells that leave thiscorridor are more likely to die, the farther they diverge from it, thehigher the rate of apoptosis gets. A family of functions is created thatdescribes this increase in mortality as cells diverge from the referencecurve

(h,η)=(c

/(η+1)³)(h−h*(η))⁴  (22)where c

is a constant, h denotes the amount of hemoglobin that a cell carries atcell age/calendric-age η and h* is the reference hemoglobin value at theage η. Polynomial functions of degree 4 are used to design thisfunction, because one of their characteristics is that they have a flatbottom (i.e. apoptosis within the corridor is low) and then increasesharply. See FIG. 16 for a graphical depiction of the function

(•,η) for a fixed calendric age η. For a graphical description of thefamily of functions

that are created, see FIG. 17. Finally, to get the rate of apoptosis αfor the precursor cells, a discrete set of N+1 equidistant data pointsμ_(i), i=0, . . . , N are defined along με└μ_(min) ^(r),μ_(max) ^(s)┘ atwhich the following calculations are done

-   -   calculate the amount of Hgb h_(i) in cells of age μ_(i) using        eq. (19), (20) and (21),    -   calculate the values of the function        at the grid points η_(i) and the corresponding Hgb value h_(i)    -   interpolate the resulting values using cubic splines.

The result is the function α(h,η) describing the rate of apoptosis ofprecursor cells shown in FIGS. 18 and 19. Note that the abovecalculations have to be repeated at every time step.

7.3. Detailed Description of the Iron Compartments.

7.3.1. Plasma (Compartment 1).

The rate at which iron enters erythroid cells at time t is given byd ₁ a ₁(t)τ(t),

where the number of transferrin receptors τ on all precursor cells is

τ(t) = ∫_(μ_(min)^(r))^(μ_(max)^(r))φ(μ^(r))r(t, μ^(r)) 𝕕μ^(r) + ∫_(μ_(min)^(s))^(μ_(max)^(s))φ(μ^(s))s(t, μ^(s)) 𝕕μ^(s),

where φ(μ) is given by eq. (17). Further, r(t,μ^(r)) and s(t,μ^(s))denote the population densities of hemoglobin synthesizing cells,erythroblasts and reticulocytes, respectively. At the same time this isthe rate at which iron leaves the plasma compartment and enters theprecursor cell compartment.

Only a very small amount of iron can be found in this pool (about 3 mg;≦1%) bound to transferrin, but it is a very important compartmentbecause of its rapid turnover rate. See Finch 1982. Iron normally turnsover at least 10 times each day. See Williams Hematology. Inflow intothe plasma compartment comes from macrophages, the storage pool, ironabsorbed through the duodenum, and via intravenously administered iron.Outflow leaves the plasma compartment to the precursor cells and thestorage compartment. The flow rates k₄₁, k_(gasto,1) and k₅₁ aredependent on the hepcidin concentration H(t) at time t. They decreasewhen the hepcidin level rises and vice versa. Inflammation and full ironstores increase the release of hepcidin from the hepatocytes in theliver, whereas a stressed erythropoiesis and depleted stores suppressthe production of hepcidin. The outflow to the precursor cellcompartment k₁₂a₁ is altered by the amount of iron needed from the cellsto synthesize hemoglobin. As erythroid cells are only able to take upiron through transferrin receptors, their need for iron can beassociated with the number of receptors expressed on precursor cells.Thus, the rate k₁₂ is given by

k₁₂ = d₁(∫_(μ_(min)^(r))^(μ_(max)^(r))φ(μ^(r))r(t, μ^(r)) 𝕕μ^(r) + ∫_(μ_(min)^(s))^(μ_(max)^(s))φ(μ^(s))s(t, μ^(s)) 𝕕μ^(s)).

For an explanation of the various terms see the statements abouterythroblasts in Section 6.2. Hence, the change of amount of iron in theplasma compartment is

${\frac{\mathbb{d}\;}{\mathbb{d}t}{a_{1}(t)}} = {{{k_{41}\left( {H(t)} \right)}{a_{4}(t)}} + {{a_{{iv},{total}}(t)}{\delta_{t_{0}}(t)}} + {{k_{{gastro},1}\left( {H(t)} \right)}{a_{gastro}(t)}} + {{k_{51}\left( {H(t)} \right)}{a_{5}(t)}} - {k_{15}{a_{1}(t)}} - {{d_{1}\left( {{\int_{\mu_{\min}^{r}}^{\mu_{\max}^{r}}{{\varphi\left( \mu^{r} \right)}{r\left( {t,\mu^{r}} \right)}\ {\mathbb{d}\mu^{r}}}} + {\int_{\mu_{\min}^{s}}^{\mu_{\max}^{s}}{{\varphi\left( \mu^{s} \right)}{s\left( {t,\mu^{s}} \right)}\ {\mathbb{d}\mu^{s}}}}} \right)}{{a_{1}(t)}.}}}$

The transfer rate k₄₁(H(t)) and k₅₁(H(t)) can be described bymonotonically decreasing sigmoidal functions

$\begin{matrix}{{{k_{41}\left( {H(t)} \right)} = {\frac{a_{6} - b_{6}}{1 + {\mathbb{e}}^{{{- k_{6}}{H{(t)}}} + c_{6}}} + b_{6}}},{and}} & (23) \\{{{k_{51}\left( {H(t)} \right)} = {\frac{a_{7} - b_{7}}{1 + {\mathbb{e}}^{{{- k_{7}}{H{(t)}}} + c_{7}}} + b_{7}}},} & (24)\end{matrix}$respectively. The parameters a₆, b₆, c₆, k₆, a₇, b₇, c₇ and k₇ arepositive constants with a₆>b₆ and a₇>b₇. The hepcidin level H(t) at timet is determined by eq. (13).7.3.2. Precursor Cells (Compartment 2).

Inflow into the precursor cell compartment comes from plasma. Oneoutflow of this compartment leaves to the erythrocytes compartment anddescribes the amount of hemoglobin carried from reticulocytes, whichenter the blood stream. Thus, the flow ƒ₂₃(•) from compartment 2 tocompartment 3, is defined by the cells that leave the bone marrow andthe amount of hemoglobin these cells are carrying. First, one has tocompute how long cells that leave the bone marrow at time t haveactually stayed there as precursor cells. Hence, one needs to solveƒ({circumflex over (t)},μ_(max) ^(s)) and one denotes the root of eq.(20) by η_(max) ^(s). Note that η_(max) ^(s) might change over time asreticulocytes are released in the blood stream faster or more slowly.k ₂₃ a ₂=ƒ₂₃(v ^(s)(E(t))s(t,μ _(max) ^(s)),h(t,η _(max) ^(s))),

where v^(s)(E(t))s(t,μ_(max) ^(s)) determines how many reticulocytes areleaving the bone marrow at a certain time t and the functionh^(s)(t,η_(max) ^(s)) indicates how much hemoglobin they carry.v^(s)(E(t)) is the maturation velocity of reticulocytes depending on EPOand s(t,μ_(max) ^(s)) is the density of the reticulocyte population withmaximal maturity μ_(max) ^(s) at time t. The function h^(s)(t,η_(max)^(s)) is given by eq. (21).

The second outflow of this compartment leaves to the macrophages and isdefined by the precursor cells that die. Hence, the flow fromcompartment 2 to compartment 4 is defined by the amount of iron thatdying cells carry, i.e.

$\begin{matrix}{{k_{24}a_{4}} = {{\int_{\mu_{\min}^{r}}^{\mu_{\max}^{r}}{{\alpha\left( {h,\mu^{r}} \right)}{h\left( {t,\mu^{r}} \right)}{r\left( {t,\mu^{r}} \right)}\ {\mathbb{d}\mu^{r}}}} + {\int_{\mu_{\min}^{s}}^{\eta_{\max}^{s}{(t)}}{{\alpha\left( {h,\eta^{s}} \right)}{h\left( {t,\eta^{s}} \right)}{s\left( {t,{\mu_{\min} + {\int_{t - \eta^{s}}^{t}{{v\left( {E(s)} \right)}\ {\mathbb{d}s}}}}} \right)}\ {{\mathbb{d}\eta^{s}}.}}}}} & (25)\end{matrix}$For the first part of the above equation, η and μ^(r) are exchangeablewithout any difficulties as, for erythroblasts, both describe the cellage of the cell, see Section 6.2.

Altogether the change of amount of iron in the iron-precursorscompartment is given by

${\frac{\mathbb{d}\;}{\mathbb{d}t}{a_{2}(t)}} = {{{d_{0}\left( {{\int_{\mu_{\min}^{r}}^{\mu_{\max}^{r}}{{\varphi\left( \mu^{r} \right)}{r\left( {t,\mu^{r}} \right)}\ {\mathbb{d}\mu^{r}}}} + {\int_{\mu_{\min}^{s}}^{\mu_{\max}^{s}}{{\varphi\left( \mu^{s} \right)}{s\left( {t,\mu^{s}} \right)}\ {\mathbb{d}\mu^{s}}}}} \right)}{a_{1}(t)}} - {f_{23}\left( {{{v^{s}\left( {E(t)} \right)}{s\left( {t,\mu_{\max}^{s}} \right)}},{h^{s}\left( {t,\mu_{\max}^{s}} \right)}} \right)} - {\int_{\mu_{\min}^{r}}^{\mu_{\max}^{r}}{{\alpha\left( {h,\mu^{r}} \right)}{h\left( {t,\mu^{r}} \right)}{r\left( {t,\mu^{r}} \right)}\ {\mathbb{d}\mu^{r}}}} - {\int_{\mu_{\min}^{s}}^{\eta_{\max}^{s}{(t)}}{{\alpha\left( {h,\eta^{s}} \right)}{h\left( {t,\eta^{s}} \right)}{s\left( {t,{\mu_{\min} + {\int_{t - \eta^{s}}^{t}{{v\left( {E(s)} \right)}\ {\mathbb{d}s}}}}} \right)}\ {{\mathbb{d}\eta^{s}}.}}}}$

Note that it is not necessary to determine the function ƒ₂₃ as, in fact,there is no need to explicitly observe the rate of change of theprecursor cell compartment. The amount of iron a₂(t) in this compartmentat time t can be directly computed from eqs. (19), (20), (21), and thepopulation equations for the erythroblasts and reticulocytes. It isdefined by the total number of precursor cells and the amount ofhemoglobin they carry:

a₂(t) = ∫_(μ_(min)^(r))^(μ_(max)^(r))α(h, μ^(r))r(t, μ^(r)) 𝕕μ^(r) + ∫_(μ_(min)^(s))^(η_(max)^(s)(t))α(h, η^(s))h(t, η^(s))s(t, μ_(min) + ∫_(t − η^(s))^(t)v(E(s)) 𝕕s) .

However, it is important to know the flux of hemoglobin to themacrophage compartment, which arises from dying cells, because this flux(see eq. (25)) is a source term in the macrophage compartment.

7.3.3. Erythrocytes (Compartment 3)

Inflow into the erythrocytes compartment comes from the precursor cellcompartment, and iron leaves the compartment to the macrophagescompartment (phagocytosis of senescent cells, random break-down of cellsand neocytolysis)

$\begin{matrix}{{{k_{34}a_{3}} = {\int_{\mu_{\min}^{m}}^{\mu_{\max}^{m}}{{h\left( {{t - \mu^{m}},\eta_{\max}} \right)}{\alpha^{m}\left( {{E(t)},\mu^{m}} \right)}{m\left( {t,\mu^{m}} \right)}\ {\mathbb{d}\mu^{m}}}}},} & (26)\end{matrix}$or is excreted via bleeding

k_(3, out)a₃ = ∫_(μ_(min)^(m))^(μ_(max)^(m))h(t − μ^(m), η_(max)^(s))α_(bleed)^(m)(t)m(t, μ^(m)) 𝕕μ^(m).

An erythrocyte of age μ^(m) left the reticulocyte population at timet−μ^(m) and, making use of Assumption 42, the amount of hemoglobin itcarries can be computed using eq. (21). Further, α^(m) denotes themortality rate of erythrocytes due to phagocytosis of senescent cells,random break-down of cells, and neocytolysis. Furthermore, h is theamount of hemoglobin carried by a cell, m(t,μ^(m)) is the populationdensity of erythrocytes with cell age μ^(m) at time t, and α_(bleed)^(m)(t) is the rate with which RBCs are lost due to bleeding.

Thus, the change of the amount of iron in the erythrocyte compartment is

${\frac{\mathbb{d}\;}{\mathbb{d}t}{a_{3}(t)}} = {{f_{23}\left( {{{v^{s}\left( {E(t)} \right)}{s\left( {t,\mu_{\max}^{s}} \right)}},{h^{s}\left( {t,\mu_{\max}^{s}} \right)}} \right)} - {\int_{\mu_{\min}^{m}}^{\mu_{\max}^{m}}{{h\left( {{t - \mu^{m}},\eta_{\max}^{s}} \right)}{\alpha^{m}\left( {{E(t)},\mu^{m}} \right)}{m\left( {t,\mu^{m}} \right)}\ {\mathbb{d}\mu^{m}}}} - {\int_{\mu_{\min}^{m}}^{\mu_{\max}^{m}}{{h\left( {{t - \mu^{m}},\eta_{\max}^{s}} \right)}{\alpha_{bleed}^{m}(t)}{m\left( {t,\mu^{m}} \right)}\ {{\mathbb{d}\mu^{m}}.}}}}$

Similar to compartment 2, it is not necessary to observe the rate ofchange of the erythrocyte cell compartment to compute the amount of irona₃(t) at time t. The quantity a₃(t) can be directly computed from thepopulation class and is defined by the total number of circulating redblood cells and the amount of hemoglobin they carry:

a₃(t) = ∫_(μ_(min)^(m))^(μ_(max)^(m))h(t − μ^(m), η_(max)^(s))m(t, μ^(m)) 𝕕μ^(m).

However, the term k₃₄a₃(t) needs to be known explicitly (see eq. (26)).This flux, which is the rate of hemoglobin released by erythrocytes,appears as a source term in the macrophage compartment.

7.3.4 Macrophages (Compartment 4)

Inflow into the macrophages compartment comes from the erythrocytescompartment (erythrocytes phagocytosed by macrophages, iron set free byrandom break-down of RBCs and shedding of hemoglobin containingvesicles) and the precursor cells compartment. Outflow of thiscompartment leaves to the plasma and to the storage compartment.Therefore, the change of iron in the compartment is

${{\frac{\mathbb{d}\;}{\mathbb{d}t}{a_{4}(t)}} = {{\int_{\mu_{\min}^{r}}^{\mu_{\max}^{r}}{{\alpha\left( {h,\mu^{r}} \right)}{h\left( {t,\mu^{r}} \right)}{r\left( {t,\mu^{r}} \right)}\ {\mathbb{d}\mu^{r}}}} + {\int_{\mu_{\min}^{s}}^{\eta_{\max}{(t)}}{{\alpha\left( {h,\eta^{s}} \right)}{h\left( {t,\eta^{s}} \right)}{s\left( {t,{\mu_{\min} + {\int_{t - \eta^{s}}^{t}{{v\left( {E(s)} \right)}\ {\mathbb{d}s}}}}} \right)}\ {\mathbb{d}\eta^{s}}}} + {\int_{\mu_{\min}^{m}}^{\mu_{\max}^{m}}{{h\left( {{t - \mu^{m}},\eta_{\max}^{s}} \right)}{\alpha^{m}\left( {{E(t)},\mu^{m}} \right)}{m\left( {t,\mu^{m}} \right)}\ {\mathbb{d}\mu^{m}}}} - {{k_{41}\left( {H(t)} \right)}{a_{4}(t)}} - {k_{45}{a_{4}(t)}}}},$where k₄₁(H(t)) is given by eq. (23).7.3.5. Storage (Compartment 5)

Inflow into the storage compartment comes from the plasma and themacrophages compartment and iron leaves this compartment to the plasmacompartment. Hence, the change of iron in the storage compartment is

${{\frac{\mathbb{d}\;}{\mathbb{d}t}{a_{5}(t)}} = {{k_{15}{a_{1}(t)}} + {k_{45}{a_{4}(t)}} - {{k_{51}\left( {H(t)} \right)}{a_{5}(t)}}}},$where k₅₁(H(t)) is given in eq. (24).7.4. Mathematical Model—Iron Model

For the convenience of the reader, all equations of the Iron Model arecollected in this section:

${{\frac{\mathbb{d}\;}{\mathbb{d}t}{a_{1}(t)}} = {{{k_{41}\left( {H(t)} \right)}{a_{4}(t)}} + {{a_{{iv},{total}}(t)}{\delta_{t_{0}}(t)}} + {{k_{{gastro},1}\left( {H(t)} \right)}{a_{gastro}(t)}} + {{k_{51}\left( {H(t)} \right)}{a_{5}(t)}} - {k_{15}{a_{1}(t)}} - {{d_{1}\left( {{\int_{\mu_{\min}^{r}}^{\mu_{\max}^{r}}{{\varphi\left( \mu^{r} \right)}{r\left( {t,\mu^{r}} \right)}\ {\mathbb{d}\mu^{r}}}} + {\int_{\mu_{\min}^{s}}^{\mu_{\max}^{s}}{{\varphi\left( \mu^{s} \right)}{s\left( {t,\mu^{s}} \right)}\ {\mathbb{d}\mu^{s}}}}} \right)}{a_{1}(t)}}}},{{a_{2}(t)} = {{\int_{\mu_{\min}^{r}}^{\mu_{\max}^{r}}{{\alpha\left( {h,\mu^{r}} \right)}{r\left( {t,\mu^{r}} \right)}\ {\mathbb{d}\mu^{r}}}} - {\int_{\mu_{\min}^{s}}^{\eta_{\max}^{s}{(t)}}{{\alpha\left( {h,\eta^{s}} \right)}{h\left( {t,\eta^{s}} \right)}{s\left( {t,{\mu_{\min} + {\int_{t - \eta^{s}}^{t}{{v\left( {E(s)} \right)}\ {\mathbb{d}s}}}}} \right)}{\mathbb{d}\eta^{s}}}}}}$  a₃(t) = ∫_(μ_(min)^(m))^(μ_(max)^(m))h(t − μ^(m), η_(max)^(s))m(t, μ^(m)) 𝕕μ^(m)${{\frac{\mathbb{d}\;}{\mathbb{d}t}{a_{4}(t)}} = {{\int_{\mu_{\min}^{r}}^{\mu_{\max}^{r}}{{\alpha\left( {h,\mu^{r}} \right)}{h\left( {t,\mu^{r}} \right)}{r\left( {t,\mu^{r}} \right)}\ {\mathbb{d}\mu^{r}}}} + {\int_{\mu_{\min}^{s}}^{\eta_{\max}^{s}{(t)}}{{\alpha\left( {h,\eta^{s}} \right)}{h\left( {t,\eta^{s}} \right)}{s\left( {t,{\mu_{\min} + {\int_{t - \eta^{s}}^{t}{{v\left( {E(s)} \right)}\ {\mathbb{d}s}}}}} \right)}\ {\mathbb{d}\eta^{s}}}} + {\int_{\mu_{\min}^{m}}^{\mu_{\max}^{m}}{{h\left( {{t - \mu^{m}},\eta_{\max}^{s}} \right)}{\alpha^{m}\left( {{E(t)},\mu^{m}} \right)}{m\left( {t,\mu^{m}} \right)}\ {\mathbb{d}\mu^{m}}}} - {{k_{41}\left( {H(t)} \right)}{a_{4}(t)}} - {k_{45}{a_{4}(t)}}}},\mspace{20mu}{{\frac{\mathbb{d}\;}{\mathbb{d}t}{a_{5}(t)}} = {{k_{15}{a_{1}(t)}} + {k_{45}{a_{4}(t)}} - {{k_{51}\left( {H(t)} \right)}{a_{5}(t)}}}},\mspace{20mu}{{\frac{\mathbb{d}\;}{\mathbb{d}t}{H(t)}} = {{\frac{1}{TBV}{H_{prod}(t)}} - {c_{\deg}^{H}{H(t)}}}},$where the initial values are given by

-   -   a₁(0)=a_(1,0),    -   a₂(0)=a_(2,0,)    -   a₃(0)=a_(3,0),    -   a₄(0)=a_(4,0),    -   a₅(0)=a_(5,0),    -   H(0)=H₀.

Moreover, a number of auxiliary equations are needed

${{k_{{gastro},1}\left( {H(t)} \right)} = {\frac{a_{5} - b_{5}}{1 + {\mathbb{e}}^{{{- k_{5}}{H{(t)}}} + c_{5}}} + b_{5}}},{{k_{41}\left( {H(t)} \right)} = {\frac{a_{6} - b_{6}}{1 + {\mathbb{e}}^{{{- k_{6}}{H{(t)}}} + c_{6}}} + b_{6}}},{{k_{51}\left( {H(t)} \right)} = {\frac{a_{7} - b_{7}}{1 + {\mathbb{e}}^{{{- k_{7}}{H{(t)}}} + c_{7}}} + b_{7}}},$where the hepcidin feedback is given by

${{{H_{prod}(t)} = {\frac{a_{H} - b_{H}}{1 + {\mathbb{e}}^{{{- k_{H}}{C{(t)}}} + c_{H}}} + b_{H}}},{{C(t)} = {{{- c_{U}}{U(t)}} + {c_{a_{5}}{a_{5}(t)}} + {c_{i}{{\mathbb{i}}(t)}}}},{{U(t)} = {{\int_{\mu_{\min}^{r}}^{\mu_{\max}^{r}}{{r\left( {t,\mu^{r}} \right)}{\mathbb{d}\mu^{r}}}} + {\int_{\mu_{\min}^{s}}^{\mu_{\max}^{s}}{{s\left( {t,\mu^{s}} \right)}\mathbb{d}}}}}}{\mu^{s}.}$

Furthermore, we have

${\varphi(\mu)} = \left\{ {{\begin{matrix}{{{3.3016\mspace{11mu}\mu^{3}} - {32.127\mspace{11mu}\mu^{2}} + {74.0952\;\mu} + 30},} & {{{{for}\mspace{14mu}\mu} \in \left\lbrack {\mu_{\min}^{r},\mu_{\max}^{r}} \right\rbrack},} \\{10,} & {{{{for}\mspace{14mu}\mu} \in \left\lbrack {\mu_{\min}^{s},\mu_{\max}^{s}} \right\rbrack},}\end{matrix}\mspace{20mu}{h\left( {t,\mu^{r}} \right)}} = {{\frac{\mathbb{d}}{TBV}{\int_{t - \mu^{r}}^{t}{{\varphi\left( {s - t + \mu^{r}} \right)}{a_{1}(s)}{\mathbb{d}s}{h\left( {t,\eta^{s}} \right)}}}} = {{h\left( {{t - \eta^{s} + \mu_{\min}^{s}},\mu_{\max}^{r}} \right)} + {\varphi\;\frac{\mathbb{d}}{TBV}{\int_{t - \eta^{s} + \mu_{\min}^{s}}^{t}{{a_{1}(s)}{{\mathbb{d}s}.}}}}}}} \right.$

The calendric age η of cells is calculated using

${{f\left( {\overset{\_}{t},\overset{\_}{\mu}} \right)} = {\mu_{\min} - \overset{\_}{\mu} + {\int_{\hat{t}}^{\overset{\_}{t}}{{v\left( {E(s)} \right)}{\mathbb{d}s}}}}},$and η=t−{circumflex over (t)}.

Further, to create the function α(h,η) (for η either μ^(r) or η^(s) canbe inserted) a discrete set of N+1 equidistant data points μ_(i), i=0, .. . , N, along με└μ_(min) ^(r),μ_(max) ^(s)┘ is defined and at everytime step the following calculations are done

-   -   calculate the amount of Hgb h_(i) in cells of age μ_(i) using        the formula for η, ƒ(t,μ), h(t,μ^(r)) and h(t,η^(s)),    -   calculate the values of the function        at the grid points η_(i) and the corresponding Hgb value h_(i)        (h,η)=(c        /(η+1)³)(h−h*(η))⁴    -   interpolate the resulting values using cubic splines.        Derivation of Population Equations        A.1. Population Equations with One Structural Variable.

Let u(t,•) denote the density of the population at time t, i.e., for0≦a<b, ∫_(a) ^(b)u(t,x)dx is the number of cells with maturity levelxε[a,b]. The rate of change is

$\begin{matrix}{{{\frac{\mathbb{d}}{\mathbb{d}t}{\int_{a}^{b}{{u\left( {t,x} \right)}{\mathbb{d}x}}}} = {\int_{a}^{b}{{u_{t}\left( {t,x} \right)}{\mathbb{d}x}}}},{t \geq 0},{0 \leq a < {b.}}} & \left( {A{.1}} \right)\end{matrix}$

Here, u is assumed to be continuously differentiable with respect to t.The rate of change with respect to t of the number of cells withmaturity level in [a,b] is the sum of several components:

-   -   the flux of cells across the boundary points a and b at time t        because of changing maturity;    -   the rate at which cells with maturity in [a,b] die at time t;    -   the rate at which cells with maturity in [a,b] proliferate.        a) Flux Across a and b.

The flux of cells across maturity level is denoted x at time t by q(t,x)and assume that the maturity of a cell is changing according to{dot over (x)}(t)=v(t,x(t)), t≧0,where the “velocity” v is continuous in t and Lipschitzean with respectto x, t≧0, 0≧0. Moreover, it is assumed that v(t,x)>0, i.e., thematurity of a cell is always strictly increasing.

Then the flux of cells across a maturity level x is given byq(t,x)=v(t,x)u(t,x), t≧0, x≧0.Consequently, the net flux of cells at time t into the interval [a,b] is

$\begin{matrix}\begin{matrix}{{{q\left( {t,a} \right)} - {q\left( {t,b} \right)}} = {{{v\left( {t,a} \right)}{u\left( {t,a} \right)}} - {{v\left( {t,b} \right)}{u\left( {t,b} \right)}}}} \\{= {- {\int_{a}^{b}\left( {{\left( {{v\left( {t,x} \right)}{u\left( {t,x} \right)}} \right)_{x}{\mathbb{d}x}},{t \geq 0.}} \right.}}}\end{matrix} & \left( {A{.2}} \right)\end{matrix}$b) Mortality.

Let x_(i) ^(N)=a+iΔx, i=0, . . . , N, with Δx=(b−a)/N be a uniform meshon [a,b]. By the mean value theorem for integrals, the number of cellswith maturity level xε[x_(i) ^(N),x_(i+1) ^(N)] at time t for somex*_(i)ε[x_(i) ^(N),x_(i+1) ^(N)] is given byu(t,x* _(i))Δx.The rate at which cells with maturity in [x_(i) ^(N),x_(i+1) ^(N)] dieat time t is assumed to be proportional to the number of cells withmaturity in this interval and therefore is given byα(t,x** _(i))u(t,x* _(i))δx,where α(t,x)≧0, t≧0, x≧0, is the mortality rate which is assumed to be acontinuous function and x**_(i)ε[x_(i) ^(N),x_(i+1) ^(N)]. Summing i=0,. . . , N−1 and taking N→∞ yields∫_(a) ^(b)α(t,x)u(t,x)dx, t≧0,  (A.3)for the rate at which cells with maturity level in [a,b] die at time t.c) Proliferation.

The assumption concerning proliferation is that a cell with maturitylevel x dividing at time t gives rise to two daughter cells withmaturity x. By analogous considerations as in the case of mortality, thefollowing expression is obtained for the rate at which cells enter theinterval [a,b] at time t because of proliferation:2∫_(a) ^(b)β(t,x)u(t,x)dx−∫ _(a) ^(b)β(t,x)u(t,x)dx=∫ _(a)^(b)β(t,x)u(t,x)dx, t≧0,  (A.4)where β(t,x)≧0 t≧0, x≧0, is the proliferation rate in the populationwhich is assumed to be continuous.d) The Model.

Using (A.1)-(A.4) yields∫_(a) ^(b) u _(t)(t,x)dx=−∫ _(a) ^(b)((v(t,x)u(t,x))_(x) dx−∫ _(a)^(b)α(t,x)u(t,x)dx+∫ _(a) ^(b)β(t,x)u(t,x)dx, t≧0, x≧0.

Assuming that the integrands in the above integrals are continuous andobserving that the interval [a,b]⊂R₊ is arbitrary yieldsu _(t)(t,x)+((v(t,x)u(t,x))_(x)=(β(t,x)−α(t,x))u(t,x), t≧0, x≧0.  (A.5)In addition to this equation, an initial condition needs to beprescribed asu(0,x)=u ₀(x), x≧0and a boundary condition for x=0. Assuming v(t,0)>0, t≧0, this boundarycondition is written as a flux condition at x=0:v(t,0)u(t,0)=ƒ₀(t), t≧0.A.2. Population Equations with Several Structural Variables.

Models for populations structured by several attributes can be derivedanalogously to those for populations with just one attribute. Forsimplicity of notation, populations are restricted to three structuralvariables and cell populations with the attributes cell age x, amount ofhemoglobin in the cell y, and amount of transferrin receptors on thecell membrane z are considered.

Let u(t,x,y,z) be the density of the cell population and Ω be a subsetof the positive cone R₊ ³ in R³ which is the admissible set for theattributes. Then the number of cells with attributes (x,y,z)εΩ at time tis given by∫_(Ω) u(t,x,y,z)dω.

Of interest is the rate at which the number of cells with attributes inΩ at time t changes, i.e., in the quantity

$\begin{matrix}{{{\frac{\mathbb{d}}{\mathbb{d}t}{\int_{\Omega}^{\;}{{u\left( {t,x,y,z} \right)}{\mathbb{d}\omega}}}} = {\int_{\Omega}^{\;}{{u_{t}\left( {t,x,y,z} \right)}{\mathbb{d}\omega}}}},} & \left( {A{.6}} \right)\end{matrix}$where it is assumed that u is continuously differentiable with respectto t. The rate at which the number of cells with attributes in Ω changesat time t is composed of several components:

-   -   the flux of cells across ∂Ω because the attributes change;    -   the rate at which cells with attributes in Ω die at time t;    -   the rate at which cells enter or leave Ω at time t because of        proliferation.        a) Flux Across ∂Ω.

The attributes of a cells change according to{dot over (x)}(t)=ƒ(t,x(t),y(t),z(t)),{dot over (y)}(t)=g(t,x(t),y(t),z(t)),{dot over (z)}(t)=h(t,x(t),y(t),z(t)),  (A.7)where ƒ, g, h are functions R₊×R₊ ³→R which are continuous with respectto t and Lipschitzean with respect to x, y, z. Then the flux of cells attime t and at (x,y,z)εR₊ ³ is given by

${F\left( {t,x,y,z} \right)} = {{u\left( {t,x,y,z} \right)}{\begin{pmatrix}{f\left( {t,x,y,z} \right)} \\{g\left( {t,x,y,z} \right)} \\{h\left( {t,x,y,z} \right)}\end{pmatrix}.}}$

Consequently the flux of cells across ∂Ω at time t is given by∫_(∂Ω)

F(t,x,y,z),v(x,y,z)

dσ,  (A.8)where v(x,y,z) denotes the unit outside normal to ∂Ω at the point(x,y,z). By Gauss' integral theorem we get

$\begin{matrix}{{\int_{\partial\Omega}{\left\langle {{F\left( {t,x,y,z} \right)},{v\left( {x,y,z} \right)}} \right\rangle{\mathbb{d}\sigma}}} = {{\int_{\Omega}{{{divF}\left( {t,x,y,z} \right)}{\mathbb{d}\omega}}} = {\int_{\Omega}{\left( {\left( {{f\left( {t,x,y,z} \right)}{u\left( {t,x,y,z} \right)}} \right)_{x} + \left( {{g\left( {t,x,y,z} \right)}{u\left( {t,x,y,z} \right)}} \right)_{y} + \left( {{h\left( {t,x,y,z} \right)}{u\left( {t,x,y,z} \right)}} \right)_{z}} \right){{\mathbb{d}\omega}.}}}}} & \left( {A{.9}} \right)\end{matrix}$b) Mortality.

Let Ω=∪_(i=1) ^(N)Ω_(i) such that diam Ω_(i)<ε, i=1, . . . , N. By themean value theorem for integrals, the number of cells with attributes inΩ_(i) at time t is given by u(t,x*_(i),y*_(i),z*_(i))ΔΩ_(i), whereΔΩ_(i) denotes the volume of Ω_(i) and (x*_(i),y*_(i),z*_(i))εΩ_(i). Itis assumed that the rate at which cells with attributes in Ω_(i) die attime t is proportional to the number of cells with attributes in Ω_(i)at time t, i.e., it is given byα(t,x** _(i) ,y** _(i) ,z** _(i))u(t,x* _(i) ,y* _(i) ,*z _(i)),where (x**_(i),y**_(i),z**_(i)), (x*_(i),y*_(i),z*_(i))εΩ_(i) andα(t,x,y,z) is the mortality rate, i.e., the rate at which cells withattributes (x,y,z) die at time t. Summing up and taking ε→0, the rate atwhich individuals with attributes in Ω die at time t is given by∫_(Ω)α(t,x,y,z)u(t,x,y,z)dω.  (A.10)c) Proliferation.

A basic assumption for proliferation is that division at time t of acell with cell age x, amount of hemoglobin y and amount of transferrinreceptors z results in two cells each with cell age x, amount ofhemoglobin y/2 and amount of transferrin receptors z/2. Furthermore, therate at which cells with attributes (x,y,z) divide at time t is assumedto be given by β(t,x,y,z).

For a subset Ω⊂R₊ ³ the subsets Ω^((1/2)), Ω⁽²⁾⊂R₊ ³ are defined byΩ^((1/2))={(x,y,z)|(x,2y,2z)εΩ}, Ω⁽²⁾={(x,y,z)|(x,y/2,z/2)εΩ}.From now on, it is assumed thatΩ^((1/2))∩Ω=Ø and Ω⁽²⁾∩Ω=Ø.  (A.11)

The following results are important when the equations governing thedevelopment of a structured population are derived:

Lemma A.1 For any (x₀,y₀,z₀)εR₊ ³ with y₀>0 or z₀>0 there exists asubset Ω⊂R₊ ³ with (x₀,y₀,z₀)εΩ satisfying condition (A.11).

Proof. The case y₀>0 is the only case considered. For ε=y₀/3>0, setΩ={(x,y,z)εR ₊ ³ ∥x−x ₀ |<ε,|y−y ₀ |<ε,|z−z ₀|<ε}and assume that (x,y,z)εΩ∩Ω^((1/2)). By the definition of Ω^((1/2)),this is equivalent to|x−x ₀ |<ε,|y−y ₀|<ε,|2y−y ₀ |<ε,|z−z ₀|<ε, and |2z−z ₀|<ε.The inequalities for y and the definition of ε imply that y>y₀−ε=2y₀/3and y<(y₀+ε)/2=2y₀/3, a contradiction, which proves that Ω^((1/2))∩Ω=Ø.The proof for Q⁽²⁾∩Ω=Ø is analogous.

The rate of change of the number of cells with attributes in Ω due toproliferation equals the rate at which the daughter cells of cells withattributes in Ω⁽²⁾ move into Ω minus the rate at which the daughtercells of cells with attributes in Ω leave Ω (and move to Ω^((1/2))). Byanalogous considerations as in the case of the effect of mortality, thefollowing expression is obtained for the rate at which the number ofcells with attributes in Ω changes because of proliferation:2∫_(Ω) ₍₂₎ β(t,{tilde over (x)},{tilde over (y)},{tilde over(z)})d{tilde over (ω)}−∫_(Ω)β(t,x,y,z)u(t,x,y,z)dω=8∫_(Ω)β(t,x,2y,2z)u(t,x,2y,2z)dω−∫_(Ω)β(t,x,y,z)u(t,x,y,z)dω  (A.12)d) The model.

Using (A.6), (A.9), (A.10) and (A.12) yields∫_(Ω) u _(t)(t,x,y,z)dω=−∫_(Ω)((ƒ(t,x,y,z)u(t,x,y,z))_(x)+(g(t,x,y,z)u(t,x,y,z))_(y)+(h(t,x,y,z)u(t,x,y,z))_(z))dω−∫_(Ω)(α(t,x,y,z)+β(t,x,y,z))u(t,x,y,z)dω+8∫_(Ω)β(t,x,2y,2z)u(t,x,2y,2z)dω.Assuming that the integrands in the above integrals are continuous andusing Lemma 1 yields

$\begin{matrix}{{{{u_{t}\left( {t,x,y,z} \right)} + \left( {{f\left( {t,x,y,z} \right)}{u\left( {t,x,y,z} \right)}} \right)_{x} + \left( {{g\left( {t,x,y,z} \right)}{u\left( {t,x,y,z} \right)}} \right)_{y} + \left( {{h\left( {t,x,y,z} \right)}{u\left( {t,x,y,z} \right)}} \right)_{z}} = {{{- \left( {{\alpha\left( {t,x,y,z} \right)} + {\beta\left( {t,x,y,z} \right)}} \right)}{u\left( {t,x,y,z} \right)}} + {8\;{\beta\left( {t,x,{2y},{2z}} \right)}{u\left( {t,x,{2y},{2z}} \right)}}}},{t \geq 0},{x > 0},{y > 0},{z > 0.}} & \left( {A{.13}} \right)\end{matrix}$An initial condition is prescribedu(0,x,y,z)=u ₀(x,y,z), x>0, y>0, z>0,and boundary conditions which, because of their special form, areconsidered directly in Section 2.A.3. List of Parameters for Models Described in Sections 6 and 7

TABLE 4 List of parameters and their meaning. Parameter Meaning β^(p)proliferation rate for BFU-E cells β^(q) proliferation rate for CFU-Ecells β^(r) proliferation rate for erythroblasts μ_(max) ^(p) maximalmaturity for BFU-E cells μ_(min) ^(q) minimal maturity for CFU-E cellsμ_(max) ^(q) maximal maturity for CFU-E cells μ_(min) ^(r) minimalmaturity for erythroblasts μ_(max) ^(r) maximal maturity forerythroblasts μ_(min) ^(s) minimal maturity for marrow reticulocytesμ_(max) ^(s) maximal maturity for marrow reticulocytes α_(rand) ^(m)intrinsic mortality rate for erythrocytes a₁, b₁, c₁, k₁ constants forthe sigmoid apoptosis rate for CFU-E cells a₂, b₂, c₂, k₂ constants forthe sigmoid maturation velocity for marrow reticulocytes d₁ rateconstant for transferrin to bind to its receptor d rate constant withwhich hemoglobin is synthesized after transferrin bound to the TfR a₃,b₃, c₃, k₃ constants for the sigmoid function governing the release ofendogenous EPO μ_(min) ^(m,n) lower bound of erythrocytes which arepossibly exposed to neocytolysis μ_(max) ^(m,n) lower bound oferythrocytes which are possibly exposed to neocytolysis μ_(max) maximallife span for erythrocytes c_(E) constant in the mortality rate forerythrocytes τ_(E) EPO threshold for neocytolysis c_(deg) ^(end)degradation rate of EPO released by the kidneys c_(deg) ^(ex)degradation rate of administered EPO (Epoetin alfa) S₀ number of cellscommitting to the erythroid lineage TBV total blood volume a₄, b₄, c₄,k₄ constants for the polynomial function describing the number of TfRexpressed on erythroblasts η calendric age of the cell η^(s) calendricage of the reticulocyte c_(x) constant in the function describing theapoptosis rate for cells with age η at time t h* reference hemoglobinvalue k_(ij) transfer rates from iron compartment i to iron compartmentj a_(iv) the amount of iron administered intravenously a_(gastro) theamount of iron in the duodenum a₅, b₅, c₅, k₅ constants for the sigmoidfunction governing the flow rate from the duodenum into the plasmacompartment a_(H), b_(H), c_(H), k_(H) constants for the sigmoidfunction governing the release of hepcidin by the liver c_(U), c_(a) ₅ ,c_(i) constants for the auxiliary equation used in the sigmoid functiongoverning the release of hepcidin by the liver c_(deg) ^(H) degradationrate of hepcidin a₆, b₆, c₆, k₆ constants for the sigmoid functiondescribing the flow rate from the macrophages into the plasmacompartment a₇, b₇, c₇, k₇ constants for the sigmoid function governingthe flow rate from the storage into the plasma compartmentA.4. Derivation of the Hgb Accumulation Function.

Let h(t,μ) respectively φ(μ) denote the amount of hemoglobin in a cellof maturity μ≧μ_(min) at time tεR respectively the number of transferrinreceptors on the membrane of a cell of maturity μ≧μ_(min). Then the rateat which h(t,μ) changes in time is given by

$\begin{matrix}{{{\frac{\mathbb{d}}{\mathbb{d}t}{h\left( {t,{\mu(t)}} \right)}} = {{{\frac{\partial}{\partial t}{h\left( {t,{\mu(t)}} \right)}} + {\frac{\partial}{\partial\mu}{h\left( {t,{\mu(t)}} \right)}{v(t)}}} = {{\mathbb{d}{\phi\left( {\mu(t)} \right)}}{a_{1}(t)}}}},\mspace{20mu}{t \in R},} & \left( {A{.4}{.1}} \right)\end{matrix}$where d>0 is a constant, a₁(t) denotes the iron concentration in plasmaat time t and v(t) is the velocity at which the maturity of a cellincreases, v(t)={dot over (μ)}(t). For equation (A.4.1) one has theboundary conditionμ(t,μ _(min))≡0,tεR.  (A.4.2)We solve the problem (A.4.1), (A.4.2) using the method ofcharacteristics. Let s→(t(s),μ(s),u(s)), s≧0, be the characteristicthrough a point (t₀,μ_(min),0). Then the functions t(s), μ(s) and u(s)are solutions of the following initial value problem (“prime(′)”=d/ds):t′(s)=1, t(0)=t ₀,μ′(s)=v(t(s)), μ(0)=μ_(min),u′(s)=dφ(μ(s))a ₁(t(s)), u(0)=0.  (A.4.3)This gives, for s≧0,t(s)=t ₀ +s,  (A.4.4)μ(s)=μ_(min)+∫₀ ^(s) v(t ₀+σ)dσ,  (A.4.5)u(s)=h(t(s), μ(s))=d∫ ₀ ^(s)φ(μ_(min)+∫₀ ^(σ) v(t ₀+ρ)dρ)a ₁(t₀+σ)dσ.  (A.4.6)

In order to determine h(t,μ) with tεR and μ≧μ_(min) one has to computetwo numbers s≧0 and {circumflex over (t)}εR such that t=t(s), μ=μ(s) andh(t,μ)=u(s) where t(•), μ(•) and u(•) solve system (A.4.3) witht₀={circumflex over (t)}. From equations (A.4.4) and (A.4.5) for s=s weget{circumflex over (t)}=t−s, μ=μ _(min)+∫₀ ^(t-{circumflex over (t)})v({circumflex over (t)}+σ)dσ=μ _(min)+∫_({circumflex over (t)}) ^(t)v(σ)dσ.This implies s=t−t and that {circumflex over (t)} has to satisfy theequationg({circumflex over (t)}):=∫_({circumflex over (t)}) ^(t) v(σ)dσ=μ−μ_(min).  (A.4.7)The function {circumflex over (t)}→g({circumflex over (t)}) ismonotonically decreasing with g(t)=0. Equation (A.4.7) has a uniquesolution {circumflex over (t)}<t if, for instance, v(t)≧δ fortε[T,{circumflex over (t)}] with δT≧μ_(max)−μ_(min) or if for some t*εRone has∫_(−∞) ^(t) ^(*) v(σ)dσ=∞.In case v≡1 one gets μ=μ_(min)+t−{circumflex over (t)} ands=t−{circumflex over (t)}, i.e.,{circumflex over (t)}=t −(μ−μ_(min)) and s=μ−μ _(min).Furthermore, one hash( t,μ )=u( s )=d∫ ₀ ^(μ−μ) ^(min) φ(μ_(min)+σ)a ₁(σ+ t−(μ−μ_(min)))dσ=d∫ _(−(μ−μ) _(min) ₎ ⁰φ(μ+σ)a ₁( t +σ)dσ.Anemia of Renal Disease

Anemia affects almost all patients with chronic kidney disease (CKD). Itis caused by failure of renal excretory and endocrine function. It oftendevelops once renal function decreases to 50% and the degree of anemiaincreases with severity of renal failure. See Strippoli, et. al. Anemiadevelops because of a deficiency in endogenous erythropoietin productionby the kidneys, increased blood losses, functional or absolute irondeficiency and decreased red cell survival. See J. B. Wish, ClinicalJournal of the American Society of Nephrology, 1:S4-S8, 2006. Concerningendogenous erythropoietin production by the kidneys, one has to notethat kidneys of end-stage renal disease patients can still produce someEPO, and thus can maintain higher hemoglobin levels than those found inanephric patients. Reasons for blood losses can be purpura,gastrointestinal and gynecologic bleeding (which occur in one third toone half of all CKD patients), frequent blood sampling, blood left inthe extracorporeal circuit, multiple vascular surgeries, etc.Furthermore, neocytolysis contributes to renal anemia and it explainsthe often demonstrable hemolytic component and the worsening ofhemolysis with more pronounced renal disease. A description ofneocytolysis is provided above. It further explains the responsivenessof hemolysis to ESA therapy and the increased efficiency of subcutaneous(s.c.) compared to intravenous therapy (i.v.), because in s.c.administration nadirs in EPO levels which precipitate neocytolysis areavoided. See Rice 1999. In general, renal anemia is normocytic andnormochromic and the number of reticulocytes is normal or slightlydecreased, which is inappropriate in the context of a reduced RBCpopulation. Certain deficiency states, especially iron, but also folateor vitamin B₁₂ deficiency may alter the nature of the anemia. See Finch1982; M. Polenakovic and A. Sikole, Journal of American Society ofNephrology, 7:1178-1182, 1996.

Untreated anemia is associated with decreased oxygen delivery to thetissues. For compensatory reasons, cardiac output increases, resultingin left ventricular hypertrophy. Cardiac disease is the most commoncause of death among patients who are on maintenance dialysis. Partialcorrection of anemia in these patients was shown to reduce cardiacischemia and ameliorate cardiomegaly, thus reducing cardiac relatedmorbidity. See A. Besarab, W. K. Bolton, J. K. Browne, J. C. Egrie, A.R. Nissenson, D. M. Okamoto, S. J. Schwab, and D. A. Goodkin, The NewEngland Journal of Medicine, 339:584-590, 1998. Further consequences ofuncorrected anemia are decreased cognition and mental acuity and overalldecrease in patient welfare.

ESA Therapy

Erythropoiesis stimulating agents have been used to treat anemia inpatients suffering from chronic renal failure for more than two decades.Optimal hemoglobin targets are still a matter of discussion. Studieshave shown (see e.g. Singh et al.; Strippoli et al., R. N. Foley, B. M.Curtis, and P. S. Parfrey. Clinical Journal of the American Society ofNephrology, 4:726-733, 2009 (hereinafter “Foley et al.”) that a partialcorrection of anemia is preferable to a full correction. A number ofcomplications can occur in patients with CKD when they havenear-normal/normal hemoglobin levels, i.e., higher vascular accessthrombosis, hypertension and greater requirements for antihypertensives,cardiovascular events, earlier need for renal replacement therapy andhigher mortality. Normal hemoglobin levels for women are in a range ofbetween about 12 g/dl and 16 g/dl, and in a range of between about 13g/dl and about 17.5 g/dl for men.

Defining an optimal hemoglobin target is not the only issue regardingESA therapy. Another problem is: how can the target hemoglobin beachieved and how to keep the patient near this hemoglobin level over along time period? The dose and frequency of administration of an ESAtreatment regimen are most often determined based on the priorexperience of the physician and on established guidelines. This approachbears some limitations and level of hemoglobin tends to fluctuategreatly and cycling phenomena are observed. An analysis of 31,267patients on hemodialysis in the Fresenius Medical Care-North Americadatabase found that only 5% of patients persistently remained within adesired Hb range of 11-12 g/dl for a period of 6 months. See A. J.Collins, R. M. Brenner, J. J. Ofman, E. M. Chi, N. Stuccio-White, M.Krishnan, C. Solid, N. J. Ofsthun, and J. M. Lazarus. American Journalof Kidney Diseases, 46:481-488, 2005 (hereinafter “Collins et al.”).Fishbane et al. analyzed data of dialysis patients collected over fiveyears in a hospital and came to a similar conclusion. See S. Fishbaneand J. S. Berns. Kidney International, 68:1337-1343, 2005 (hereinafter“Fishbane et al.”). More than 90% of the patients experienced hemoglobincycling. The authors state that changes in ESA dose were the mostimportant driver and were associated with hemoglobin excursion in about80% of cases. The ESA dose adjustment protocol that was used was similarto the protocol used in most dialysis centers.

Therapy with ESAs is quite different than biological erythropoietinsecretion. In hemodialysis patients, i.v. administration route isprimarily used, because of the availability of venous access.Intravenous treatment involves short, intermittent, non-physiologicbursts of EPO concentration. The bursts are followed by a fast declineto very low levels of EPO. These fluctuations in plasma concentration donot coincide, either temporally or in magnitude, with physiologicperturbations. Therefore, it may not be surprising that Hb levelsfluctuate widely and that it is extremely difficult for physicians toadjust dosing schemes such that no cycling phenomena occur. Note thatfluctuations in hemoglobin result in varying oxygen delivery to vitalorgans. Consequences include repeated episodes of relative ischemia andcompensatory mechanisms in organs (e.g., heart) that may result indisordered growth signals, pathologic organ function and worsenedpatient outcomes. (See Fishbane et al.).

A predictive model of erythropoiesis can help deal with this situation.For instance, a physician can try different ESA treatment regimens andobserve their effects on hemoglobin levels over the next few months.Dosing regimens can be tested/chosen with regard to avoidance ofneocytolysis, minimal amounts of ESA administered and avoidance ofcycling patterns in Hb concentration.

Iron Therapy

Patients with anemia of chronic kidney disease have to be followed forsymptoms of iron deficiency. See J. W. Fisher. Experimental Biology andMedicine, 28:1-24, 2003. 80-90% of dialysis patients on ESA therapy willrequire iron at some stage. See R. M. Schaefer and L. Schaefer.Nephrology Dialysis Transplantation, 13:9-12, 1998. This very pronouncedneed for supplementary iron has different reasons. On one hand ironstores can be depleted (absolute iron deficiency). Iron loss inhemodialysis patients (due to continuous gastrointestinal, purpural andgynecological bleeding, frequent blood sampling, surgeries, . . . ) isabout 1500-3000 mg/year, as compared to iron loss for a healthy adultthat is about 400-800 mg/year, and can be even higher under certaincircumstances. Hence, the daily need for iron can be well above theabsorptive capacity of the intestinal. This is aggravated by the factthat the uptake via the duodenum is often impaired in these patients. Onthe other hand, a functional iron deficiency is often observed in renalanemia. During ESA therapy, the number of erythroid progenitor andprecursor cells in the bone marrow increase drastically and this imposesa lot of stress on systemic iron homeostasis. Supply and demand often donot match. It is difficult, even for healthy persons, to increase therate at which iron is released from stores and is recycled fromhemoglobin to deliver enough iron to the bone marrow to keep up withsupraphysiologic rates of RBC production during ESA treatment. Mattersare further complicated in chronic kidney disease because ofinflammation, which frequently occurs with various degrees of severity.Thus, iron utilization is regularly decreased in renal insufficiency.See above for a detailed description of the effects of inflammation onsystemic iron homeostasis.

In the pre-ESA era iron overload, because of excessive bloodtransfusions, was a major cause of morbidity in dialysis patients. Itssignificance in the post-EPO era remains unclear. See K. Kalantar-Zadeh,D. L. Regidor, C. J. McAllister, B. Michael, and Warnock D. G., Journalof American Society of Nephrology, 16:3070-3080, 2005. A continuousadministration of i.v. iron certainly results in overfilled iron-storesand imposes serious health risks. Hence, decisions on when to supply thepatient with i.v. iron and when to withdraw therapy have to bethoroughly evaluated. A mathematical model can help keep track of thecurrent iron status of a patient and can help to make decisions onadaptations of treatment.

Finally, there is a very complex interaction between hemoglobin cyclingand iron storage whose dynamical behavior under ESA and iron treatment,even for an experienced physician, is barely predictable. In Alfrey etal., the authors suggest that: “[ . . . ] the current therapeuticparadigm of hemoglobin monitoring, iron treatment, and rHuEPO treatmentresults in recurrent nonphysiologic cycling of hemoglobin levels inhemodialysis patients.”

Computation of Hematocrit and Hemoglobin Concentrations

In order to compute the hematocrit (HCT) and the hemoglobinconcentration (Hb) for a subject from the model output (which is thenumber of red blood cells circulating in blood), estimates of the totalblood volume of the subject are needed. The Nadler and Allen formula(see Nadler, S. B., Hidalgo, J. U., and Bloch, T., Surgery, 1962, Vol.51, pp. 224-232) can be used to estimate the total blood volume (TBV) ofa healthy subject according to her/his weight and height:Women: TBV [ml]=183.3+356.1×(height [m])³+33.08×weight [kg],Men: TBV [ml]=604.1+366.9×(height [m])³+32.19×weight [kg].

The TBV for a dialysis patient can be measured by radio-labeling redblood cells with chromium-51 as described above. Using these estimatesfor the total blood volume, the hematocrit (HCT) and the hemoglobinconcentration (Hb) can be computed for a patient from the number of redblood cells circulating in blood via the following formulae:

${{{HCT}\lbrack\%\rbrack} = \frac{\left( {{M(t)} \times {{MCV}\lbrack{fl}\rbrack}} \right)}{{TBV}\lbrack{ml}\rbrack}},$where M(t)=∫₀ ^(μ) ^(max) ^(m) m(t,μ^(m))dμ^(m) is the total number oferythrocytes circulating in blood and MCV is the mean corpuscular volumeof a RBC which is obtained from measurement, and

${{{Hb}\left\lbrack \frac{g}{l} \right\rbrack} = {1000 \times \frac{M(t) \times {{MCH}\lbrack{pg}\rbrack}}{{TBV}\lbrack{ml}\rbrack}}},$where MCH is the mean cellular hemoglobin, which is also obtained viameasurements.Implementation in a Computer Network

FIG. 20 illustrates a computer network or similar digital processingenvironment in which the present invention can be implemented.

Computer(s)/devices 50 and server computer(s) 60 provide processing,storage, and input/output devices executing application programs and thelike. Computer(s)/devices 50 can also be linked through communicationsnetwork 70 to other computing devices, including other devices/processes50, digital processor dialysis machines 50A, including machines withintegrated modular drug delivery devices (as shown in FIGS. 2 and 3),and server computer(s) 60. Communications network 70 can be part of aremote access network, a global network (e.g., the Internet), aworldwide collection of computers, Local area or Wide area networks, andgateways that currently use respective protocols (TCP/IP, Bluetooth,etc.) to communicate with one another. Other electronic device/computernetwork architectures are suitable.

FIG. 21 is a diagram of the internal structure of a computer (e.g.,processor/device 50, digital processor dialysis machines 50A, or servercomputers 60) in the computer system of FIG. 20. Each computer 50, 50A,60 contains a system bus 79, where a bus is a set of hardware lines usedfor data transfer among the components of a computer or processingsystem. Bus 79 is essentially a shared conduit that connects differentelements of a computer system (e.g., processor, disk storage, memory,input/output ports, network ports, etc.) that enables the transfer ofinformation between the elements which could be data collected during atreatment, but also instructions to a machine with integrated modulardrug delivery devices to follow a drug delivery scheme provided by acentral erythropoiesis modeling engine. Attached to system bus 79 is I/Odevice interface 82 for connecting various input and output devices(e.g., keyboard, mouse, displays, printers, speakers, etc.) to thecomputer 50, 60. Network interface 86 allows the computer to connect tovarious other devices attached to a network (e.g., network 70 of FIG.20). Memory 90 provides volatile storage for computer softwareinstructions 92 and data 94 used to implement an embodiment of thepresent invention (e.g., Eqs. described above and listed in Section 4 orany other erythropoiesis modeling engine detailed above). Disk storage95 provides non-volatile storage for computer software instructions 92and data 94 used to implement an embodiment of the present invention.Central processor unit 84 is also attached to system bus 79 and providesfor the execution of computer instructions.

In one embodiment, the processor routines 92 and data 94 are a computerprogram product (generally referenced 92), including a computer readablemedium (e.g., a removable storage medium such as one or more DVD-ROM's,CD-ROM's, diskettes, tapes, etc.) that provides at least a portion ofthe software instructions for the system. Computer program product 92can be installed by any suitable software installation procedure, as iswell known in the art. In another embodiment, at least a portion of thesoftware instructions can also be downloaded over a cable, communicationand/or wireless connection. In other embodiments, the programs are acomputer program propagated signal product 107 embodied on a propagatedsignal on a propagation medium (e.g., a radio wave, an infrared wave, alaser wave, a sound wave, or an electrical wave propagated over a globalnetwork such as the Internet, or other network(s)). Such carrier mediumor signals provide at least a portion of the software instructions forthe present routines/program 92.

In alternate embodiments, the propagated signal is an analog carrierwave or digital signal carried on the propagated medium. For example,the propagated signal may be a digitized signal propagated over a globalnetwork (e.g., the Internet), a telecommunications network, or othernetwork. In one embodiment, the propagated signal is a signal that istransmitted over the propagation medium over a period of time, such asthe instructions for a software application sent in packets over anetwork over a period of milliseconds, seconds, minutes, or longer. Inanother embodiment, the computer readable medium of computer programproduct 92 is a propagation medium that the computer system 50 canreceive and read, such as by receiving the propagation medium andidentifying a propagated signal embodied in the propagation medium, asdescribed above for computer program propagated signal product.

Generally speaking, the term “carrier medium” or transient carrierencompasses the foregoing transient signals, propagated signals,propagated medium, storage medium and the like.

The relevant teachings of all patents, published patent applications andliterature references cited herein are incorporated by reference intheir entirety.

While this invention has been particularly shown and described withreferences to example embodiments thereof, it will be understood bythose skilled in the art that various changes in form and details may bemade therein without departing from the scope of the inventionencompassed by the appended claims.

What is claimed is:
 1. A hemodialysis system comprising: a hemodialysismachine comprising a blood pump; a modular drug delivery devicecomprising at least one erythropoiesis stimulating agent (ESA) pump, atleast one iron replacement product pump, at least one ESA vial holder,and at least one iron replacement product vial holder; and a housing towhich the blood pump and the modular drug delivery device are secured; ablood line set comprising multiple blood lines and a drip chamber influid communication with the multiple blood lines, at least one of theblood lines being connected to the drip chamber, the blood line setbeing connected to the blood pump in a manner such that, when the bloodline set is connected to a patient and the blood pump is operated, bloodof the patient is passed through the blood line set; a drugadministration fluid line set comprising multiple drug lines, at leasttwo of the drug lines of the drug line set being connected to the dripchamber of the blood line set, and the drug line set being connected tothe at least one ESA pump and the at least one iron replacement productpump in a manner such that, when the drug line set is fluidly connectedto an ESA vial contained in the at least one ESA vial holder and to aniron replacement product vial contained in the at least one ironreplacement product vial holder, and the at least one ESA pump and theat least one replacement product pump are operated, ESA and ironreplacement product are delivered from the ESA vial and the ironreplacement vial, respectively, to the drip chamber of the blood lineset via the drug line set; and a control unit comprising a mathematicalmodel of erythropoiesis and iron metabolism and instructions that, whenexecuted, cause the control unit to adjust the patient's undesiredhematocrit and/or hemoglobin concentration to a value within a desiredrange at a predetermined time by i) employing the model to predict thepatient's hematocrit and/or hemoglobin concentration at a predeterminedtime, the model including iron homeostasis, and by ii) operating theblood pump, the at least one ESA pump, and the at least one ironreplacement product pump to deliver blood, ESA, and iron replacementproduct to the drip chamber via the at least one of the blood linesconnected to the drip chamber and the at least two of the drug linesconnected to the drip chamber, according to an ESA administrationregimen predicted to adjust the patient's hematocrit and/or hemoglobinconcentration to a value within the desired range at the predeterminedtime, wherein the mathematical model of erythropoiesis and ironmetabolism comprises the functions:$\mspace{20mu}{{{{\frac{\partial}{\partial t}{p\left( {t,\mu^{p}} \right)}} + {\frac{\partial}{\partial\mu^{p}}{p\left( {t,\mu^{p}} \right)}}} = {\beta^{p}{p\left( {t,\mu^{p}} \right)}}},\mspace{20mu}{{{\frac{\partial}{\partial t}{q\left( {t,\mu^{q}} \right)}} + {\frac{\partial}{\partial\mu^{q}}{q\left( {t,\mu^{q}} \right)}}} = {\left( {\beta^{q} - {\alpha^{q}\left( {E(t)} \right)}} \right){q\left( {t,\mu^{q}} \right)}}},\mspace{20mu}{{{\frac{\partial}{\partial t}{r\left( {t,\mu^{r}} \right)}} + {\frac{\partial}{\partial\mu^{r}}{r\left( {t,\mu^{r}} \right)}}} = {\beta^{r}{r\left( {t,\mu^{r}} \right)}}},\mspace{20mu}{{{\frac{\partial}{\partial t}{s\left( {t,\mu^{s}} \right)}} + {{v^{s}\left( {E(t)} \right)}\frac{\partial}{\partial\mu^{s}}{s\left( {t,\mu^{s}} \right)}}} = {{- \alpha^{s}}{s\left( {t,\mu^{s}} \right)}}},\mspace{20mu}{{{\frac{\partial}{\partial t}{m\left( {t,\mu^{m}} \right)}} + {\frac{\partial}{\partial\mu^{m}}{m\left( {t,\mu^{m}} \right)}}} = {{- {\alpha^{m}\left( {{E(t)},\mu^{m}} \right)}}{m\left( {t,\mu^{m}} \right)}}},\mspace{20mu}{{\frac{\mathbb{d}}{\mathbb{d}t}{E^{end}(t)}} = {{\frac{1}{TBV}{E_{in}^{end}(t)}} - {c_{\deg}^{end}{E^{end}(t)}}}},\mspace{20mu}{{\frac{\mathbb{d}}{\mathbb{d}t}{E^{ex}(t)}} = {{\frac{1}{TBV}{E_{in}^{ex}(t)}} - {c_{\deg}^{ex}{E^{ex}(t)}}}},\mspace{20mu}{where}}$$\mspace{20mu}{{{\alpha^{q}\left( {E(t)} \right)} = {\frac{a_{1} - b_{1}}{1 + {\mathbb{e}}^{{k_{1}{E{(t)}}} - c_{1}}} + b_{1}}},\mspace{20mu}{{v^{s}\left( {E(t)} \right)} = {\frac{a_{2} - b_{2}}{1 + {\mathbb{e}}^{{{- k_{2}}{E{(t)}}} + c_{2}}} + {b_{2}.\mspace{20mu}{and}}}}}$${\alpha^{m}\left( {{E(t)},\mu^{m}} \right)} = \left\{ {{\begin{matrix}{{{\gamma^{m}\left( \mu^{m} \right)} + {\min\left( {\frac{c_{E}}{{E(t)}^{k_{E}}},b_{E}} \right)}},} & {{{{forE}(t)} < \tau_{E}},{\mu_{\min}^{m,n} \leq \mu^{m} \leq \mu_{\max}^{m,n}},} \\{{\gamma^{m}\left( \mu^{m} \right)},} & {{otherwise},}\end{matrix}\mspace{20mu}{with}{\gamma^{m}\left( \mu^{m} \right)}} = \left\{ {{{\begin{matrix}\alpha_{rand}^{m} & {{{{for}\;\mu^{m}} \in \left\lbrack {\mu_{\min}^{m},{\mu_{\max}^{m} - \delta}} \right\rbrack},} \\\frac{3\;\alpha_{rand}^{m}\delta^{2}}{\left( \mu^{m} \right)^{2} - {2\left( {\mu_{\max}^{m} + \delta} \right)\mu^{m}} + {\left( {\mu_{\max}^{m} + {2\;\delta}} \right)\mu_{\max}^{m}}} & {{{{for}\;\mu^{m}} \in \left\lbrack {{\mu_{\max}^{m} - \delta},\mu_{\max}^{m}} \right)},} \\\infty & {{{for}\;\mu^{m}} \geq {\mu_{\max}^{m}.}}\end{matrix}\mspace{20mu}{and}\begin{matrix}{{\frac{\mathbb{d}}{\mathbb{d}t}{a_{1}(t)}} = {{{k_{41}\left( {H(t)} \right)}{a_{4}(t)}} + {{a_{{iv},{total}}(t)}{\delta_{t_{0}}(t)}} + {{k_{{gastro},1}\left( {H(t)} \right)}{a_{gastro}(t)}} +}} \\{{{k_{51}\left( {H(t)} \right)}{a_{5}(t)}} - {k_{15}a_{1}(t)} -} \\{{{d_{1}\left( {{\int_{\mu_{\min}^{r}}^{\mu_{\max}^{r}}{{\varphi\left( \mu^{r} \right)}{r\left( {t,\mu^{r}} \right)}{\mathbb{d}\mu^{r}}}} + {\int_{\mu_{\min}^{s}}^{\mu_{\max}^{s}}{{\varphi\left( \mu^{s} \right)}{s\left( {t,\mu^{s}} \right)}{\mathbb{d}\mu^{s}}}}} \right)}{a_{1}(t)}},}\end{matrix}{a_{2}(t)}} = {{{\int_{\mu_{\min}^{r}}^{\mu_{\max}^{r}}{{\alpha\left( {h,\mu^{r}} \right)}{r\left( {t,\mu^{r}} \right)}{\mathbb{d}\mu^{r}}}} + {\int_{\mu_{\min}^{s}}^{\eta_{\max}^{s}{(t)}}{{\alpha\left( {h,\eta^{s}} \right)}{h\left( {t,\eta^{s}} \right)}{s\left( {t,{\mu_{\min} + {\int_{t - \eta^{s}}^{t}{{v\left( {E(s)} \right)}{\mathbb{d}s}}}}} \right)}{\mathbb{d}\eta^{s}}\mspace{20mu}{a_{3}(t)}}}} = {{\int_{\mu_{\min}^{m}}^{\mu_{\max}^{m}}{{h\left( {{t - \mu^{m}},\eta_{\max}^{s}} \right)}{m\left( {t,\mu^{m}} \right)}{\mathbb{d}\mu^{m}}\begin{matrix}{{\frac{\mathbb{d}}{\mathbb{d}t}{a_{4}(t)}} = {{\int_{\mu_{\min}^{r}}^{\mu_{\max}^{r}}{{\alpha\left( {h,\mu^{r}} \right)}{h\left( {t,\mu^{r}} \right)}{r\left( {t,\mu^{r}} \right)}{\mathbb{d}\mu^{r}}}} +}} \\{{\int_{\mu_{\min}^{s}}^{\eta_{\max}^{s}{(t)}}{{\alpha\left( {h,\eta^{s}} \right)}{h\left( {t,\eta^{s}} \right)}{s\left( {t,{\mu_{\min} + {\int_{t - \eta^{s}}^{t}{{v\left( {E(s)} \right)}{\mathbb{d}s}}}}} \right)}{\mathbb{d}\eta^{s}}}} +} \\{{\int_{\mu_{\min}^{m}}^{\mu_{\max}^{m}}{{h\left( {{t - \mu^{m}},\eta_{\max}^{s}} \right)}{\alpha^{m}\left( {{E(t)},\mu^{m}} \right)}{m\left( {t,\mu^{m}} \right)}{\mathbb{d}\mu^{m}}}} -} \\{{{{k_{41}\left( {H(t)} \right)}{a_{4}(t)}} - {k_{45}{a_{4}(t)}}},}\end{matrix}\mspace{20mu}\frac{\mathbb{d}}{\mathbb{d}t}{a_{5}(t)}}} = {{k_{15}{a_{1}(t)}} + {k_{45}{a_{4}(t)}} - {{k_{51}\left( {H(t)} \right)}{a_{5}(t)}}}}}},\mspace{20mu}{{{\frac{\mathbb{d}}{\mathbb{d}t}{H(t)}} = {{\frac{1}{TBV}{H_{prod}(t)}} - {c_{\deg}^{H}{H(t)}}}};}} \right.} \right.$wherein TBV is a total blood volume, p(t,μ^(p)) is a population densityof a BFU-E cell class at time t with maturity μ^(p), β^(p) is aproliferation rate for BFU-E cells, q(t,μ^(q)) is a population densityof a CFU-E cell class at time t with maturity μ^(q), β^(q) is aproliferation rate for CFU-E cells, α^(q)(E(t)) is an apoptosis rate forCFU-E cells which depends on erythropoietin (EPO)-concentration E(t) attime t, r(t,μ^(r)) is a population density of an erythroblast class attime t with maturity μ^(r), β^(r) is a proliferation rate forerythroblasts, s(t,μ^(s)) is a population density of a marrowreticulocytes class at time t with maturity μ^(s), v^(s)(E(t)) is amaturation velocity of cells leaving the erythroblast class an enteringthe reticulocytes population class which depends on EPO-concentrationE(t) at time t, α^(s) is an apoptosis rate for marrow reticulocytes,m(t,μ^(m)) is a population density of an erythrocyte class at time twith maturity μ^(m), α^(m)(E(t),μ^(m)) is an apoptosis rate forerythrocytes which depends on EPO-concentration E(t) at time t andmaturity μ^(m), E^(end)(t) is an endogenous EPO concentration at time t,E_(in) ^(end)(t) is an amount of EPO released by a patient's kidneys,c_(deg) ^(end) is a degradation rate of endogenous EPO, E^(ex)(t) is anexogenous EPO concentration at time t, E_(in) ^(ex)(t) is an amount ofEPO administered to the patient, c_(deg) ^(ex) is a degradation rate ofexogenous EPO, a₁, b₁, c₁ and k₁ are constants for the apoptosis ratefor CFU-E cells, a₂, b₂, c₂ and k₂ are constants for the maturationvelocity for marrow reticulocytes, γ^(m)(μ^(m)) is a mortality rate oferythrocytes which depends on maturity μ^(m), b_(E), c_(E), and k_(E)are constants in the mortality rate of erythrocytes, τ_(E) is athreshold beneath which neocytolysis is triggered, μ_(min) ^(m,n) is alower bound of erythrocytes which are possibly exposed to neocytolysis,μ_(max) ^(m,n) is an upper bound of erythrocytes which are possiblyexposed to neocytolysis, α_(rand) ^(m) is an intrinsic mortality ratefor erythrocytes, δ_(t0)(t) is a Dirac delta impulse at a time ofadministration t₀, H(t) is a hepcidin concentration and time t,α_(iv,total)(t) is a total amount of iron intravenously administered attime t, α_(gastro)(t) is an amount of iron in the duodenum of thepatient at time t, k₁₅ is a transfer rate of iron from plasma tostorage, k₄₁ is a transfer rate of iron from macrophages to plasma, k₄₅is a transfer rate of iron from macrophages to plasma, k₅₁ is a transferrate of iron from storage to plasma, k_(gastro) is a transfer rate ofiron in the duodenum of the patient, a₁(t) is a concentration of iron inplasma at time t, a₂(t) is a concentration of iron in a precursor cellcompartment at time t, a₃(t) is a concentration of iron in a erythrocytecell compartment at time t, a₄(t) is a concentration of iron in amacrophage compartment at time t, a₅(t) is a concentration of iron in astorage compartment at time t, φ(μ) is a function describing a number ofTfR on a cell with cell age μ, r(t,μ^(r)) is the density oferythroblasts, s(t,μ^(s)) is the density of bone marrow reticulocytes,h(t,μ(t)) denotes an amount of hemoglobin in a cell with age μ(t) attime t, and α(h,η) is a rate of apoptosis of precursor cells withcalendric age η.
 2. The hemodialysis system of claim 1, wherein the ESAadministration regimen includes iron administration, and the modelincludes iron homeostasis with iron administration.
 3. The hemodialysissystem of claim 1, wherein predicting the patient's hematocrit and/orhemoglobin concentration includes modeling influence of iron homeostasison erythropoiesis.
 4. The hemodialysis system of claim 1, wherein thepredetermined time is in a range of between about 5 days and about 200days into the ESA administration regimen.
 5. The hemodialysis system ofclaim 1, wherein the patient undergoes a medical procedure prior,during, or after initiation of the ESA administration regimen, themedical procedure selected from the group consisting of blood donation,surgery, and dialysis, or any combination thereof.
 6. The hemodialysissystem of claim 1, wherein the patient is a dialysis patient and thedesired hematocrit is in the range of between about 28 percent and about36 percent and the desired hemoglobin concentration is in a range ofbetween about 9.5 g/dL and about 12 g/dL.